﻿Ion Necoara 2014 Personal Curs: Prof Dr Ion Necoara (ion necoara@acse pub ro) Seminar: Drd Andrei Patrascu (andrei patrascu@acse pub ro) Examinare 15 puncte seminar, lab (prezenta obligatorie) si curs 25 puncte colocviu (p^ saptamana 8) 60 puncte examen final In 1744 Leonhard Euler, unul din cei mai mari matematicieni: “Namely, because the shape of the whole universe is most perfect and, in fact, designed by the wisest creator, nothing in all of the world will occur in which no maximum or minimum rule is somehow shining forth " In latina: “ nihi omnino in mundo contingint, in quo non maximi minimive ratio quapiam eluceat " De fapt, conceptul de derivata a functiei (una din cele mai importante constructii matematice) a fost introdus de Fermat in 1637 cu scopul rezolvarii unei probleme de optimizare Note de curs/lab: Carte: Necoara, Metode de Optimizare Numerica, Politehnica Press Culegere: Necoara, Metode de Optimizare Numerica: Culegere de Probleme, Politehnica Press Alte cârti: Bertsekas, Nonlinear Programming, Athena Scientific Nesterov, Introductory Lectures on Convex Optimization: A Basic Course, Kluwer Luenberger, Linear and Nonlinear Programming, Kluwer Obiective recunoasterea/formularea de probleme de optimizare dezvoltarea de cod Matlab pentru rezolvarea problemelor de optimizare caracterizarea/identificarea soluțiilor caracterizarea limitelor de performanta a algoritmilor numerici de optimizare Subiecte atinse mulțimi si funcții convexe clase de probleme de optimizare condiții de optimalitate algorimi numerici de optimizare si proprietățile lor exemple si aplicații In cadrul cursului vom utiliza următoarele notatii: Vectori (considerați intotdeauna vector coloana) cu litere mici, xi i e x G Rn, x = : Xn Produs scalar in spațiul Euclidian: (x,y) = xTy = 52"=1x,y/ Norma Euclidiană standard ||x|| = a/(x, x) = xl + ■ ■ ■ + x^ Mulțimi cu litere mari: S, Q, U C Rn ( - orthantul nenegativ, Sț - mulțimea matricelor pozitiv semidefinite) Matrice cu litere mari: A, B, C, H E Rnxm Norma spectrala a unei matrici ||4|| = y/Xmax(ATA') Matrice pozitiv definita: A >- 0, si pozitiv semidefinita /IM In cadrul cursului vom utiliza următoarele notatii: Funcții cu litere mici: f, g, h: Rn —> R Gradientul, respectiv matricea Hessiana a funcției continuu diferentiabile f: Rn -G R: Vf(x) G Rn, V2f(x) G Rnxn Hessiana este matrice simetrica (i e in Sn)l dxi 92xi Vf(x) = df(x) 9xn _ , V2f(x) = d2f(x) _dxndxi 92fM' dxidxn 92f(x) d2x„ - Teorema de medie: fie funcția g : R —> R atunci exista 0 G [a, b] a i g(b) - g(a) = g'^b- a) = // g'(r)dT pentru f : Rn —> R luam g(t) = f(x + t(y — x)): f (y) = f (x) + (Vf(x + 6*(y - x)), y - x) 0 G /r(y) = /rW+/ (Vf(x + r(y-x)),y-x)dT Optimizarea matematica = = Urcarea unui munte f* = min f(x) xeR" s l : gi(x) R constrângeri de inegalitate si egalitate g,-, hj: Rn —> R In soluția optima x* funcția obiectiv ia valoarea cea mai mica (valoarea optima f*) in raport cu toti vectorii ce satisfac constrângerile (restricțiile) Scurta istorie: Antichitate: abordare geometrica (aplicații izoperimetrice -Problema Dido sec IX i H ) Evul mediu: abordare algebrica (folosind e g inegalitatea mediilor) După sec XIX: abordare pe baza calculului diferențial teorie: condiții de optimalitate, analiza convexa (1900-1970) algoritmi: 1947: algoritmul simplex pentru programare liniara (Dantzig) 1950: programarea dinamica (Bellman ’52) 1960: metoda gradientilor conjugați (Fletcher & Reeves ’64); metode quasi-Newton (David-Fletcher-Powel ’59) 1970: metode de penalitate si bariera (Fiacco & McCormick); programare patratica secvențiala; probleme de mari dimensiuni (Lasdon ’70); metoda subgradient (școala ruseasca) 1980-2000: metode de punct interior cu complexitate polinomiala (Karmakar '84; Nesterov & Nemirovski ’94) aplicații: inainte de 1990: cercetări operaționale (finanțe, logistica, manufacturing, scheduling) după 1990: aplicații inginerești (control, procesare de semnal, motoare de cautare, rețele); medicina (tomografie); biologie (recunoașterea de gene); fizica (predictia vremii) Problema de optimizare neliniara: f* = min f(x) xgR" s l : gi(x) R si variabila de decizie xGR" Notatii compacte: g : R" —> Rm, g(x) = [gi(x) gm(x)]r /?:R''->RP r forma compacta a problemei de optimizare neliniara ((NLP) -nonlinear programming) (NLP) : min f(x) xeR" s L: g(x) 0, xf + xj — 1 R3, unde gi(x) = -Xi, g2(x) = -X2 si g3(x) = xf +x2 - 1 o singură constrângere de egalitate definita de funcția /)(x) = X1X2 — 1 Alte definiții esențiale: Mulțimea {xgR": f(x) = c} este mulțimea nivel (conturul) a funcției f pentru valoarea c G R Mulțimea fezabila a problemei de optimizare (NLP) este: X = {x G Rn : g(x) R este continua atunci exista un punct de minim global pentru problema de minimizare minxex f(x) Exemplu: pentru funcția f(x) din figura consideram mulțimea fezabila compacta X = {x : a R, f(x) = x® + x24 + x3 g : R3 —> R, g(x) = 5xi + 6x2 h : R3 R2, /?(x) = *1 + *3 - 4 xi + x2 + 2 Daca f(x),g(x) si h(x) din problema NLP sunt afine (f(x) = cTx, g(x) = Cx — d, h(x) = Ax — b), atunci problema NLP devine o problema de programare liniara (LP sau linear programming): (LP) : min cTx xgR" s L: Cx — d 0, X3 > 0 *1 + X2 = 2 x2 + x3 = 1 xi Observam: f(x) = 3xi + 2x2 + 4x3 = X2 Constrângeri de egalitate: /i(x) : R3 —> R2 /i(x) = 4x - b = 0 o 1 1 1 0 1 *1 xj 2 1 = 0 A b Astfel, problema (1) poate fi scrisa ca problema LP standard Daca g(x) si /)(x) din problema NLP sunt afine iar funcția obiectiv f(x) este o funcție patratica, atunci problema NLP devine o problema de programare patratica (QP sau quadratic programming): (QP): 1 T min —x Qx + q x + r xeR" 2 s l : Cx — d 0 QP-urile se rezolva cu metode punct interior, “active set” sau de tip gradient: astazi se pot rezolva probleme de dimensiuni foarte mari (109 variabile) ce apar in multe aplicații: optimizarea motoarelor de cautare, procesarea de imagine, predictie meteo ( ) Exemple de probleme QP: min x? + x/ + X1X2 + 3xi + 2x2 xeR" s l : x > 0, xi + X2 = 5 Observam f (x) = xf + x/ + X1X2 + 3xi + 2x2 *1 X2 = -2 H Constrângeri de inegalitate: —l?x - 0 iar f(x) este convexa • Ierarhizarea in motoarele de cautare se pune ca o problema QP min ||Ex — x||2 xeR": eTx=l, x>0 E matrice de adiacenta (http://acse pub ro/person/ion-necoara) Probleme de optimizare convexa: (CP) : min f(x) xeR" s L: g(x) 0) sunt probleme de optimizare convexa Problemele convexe sunt, in general, mai ușor de rezolvat decât (NLP)-urile Foarte multe aplicații: procesare de semnal, control optimal, design optim • Observație: o problema de optimizare convexa unde f(x) = ^xtQx + qTx si g/(x) = + q[x + q, i e f si funcțiile din constrângerile de inegalitate sunt patratice, se numește QCQP(quadratically constrained quadratic program) Problema SDP (semidefinite programming): probleme convexe unde mulțimea fezabila este descrisa de LMI-uri Aceste probleme se numesc SDP-uri deoarece se impune ca anumite matrice sa ramana pozitiv definite O problema generala SDP poate fi formulata drept: (SDP) : min cTx xCR" n s L: Ao + A/x/=^0, Ax — b = 0, i=i unde matricele A; G Sm oricare ar fi i = 0, , n Remarcam ca problemele LP, QP, si QCQP pot fi de asemenea formulate ca probleme SDP Foarte multe aplicații: control optimal si robust, stabilitatea sistemelor dinamice, probleme de localizare, matrix completion, Problemele convexe sunt mult mai ușor de rezolvat decât cele neconvexe! (software gratuit de optimizare bun pentru probleme convexe este ) Exemplu de problema QCQP convexa scrisa ca un SDP: 1 t T min —x Qx + q x xgR" 2 s l : ~xTQix+qfx+ri ±xtQx + q |xrQix+q1r x+ri SDP echivalent min t xeR",teR In x 1/V2 = 0 l/V2xTQl l/V2Qlx -qȚx- ri >= 0 APLICAȚII REALE • se dau funcțiile de baza fi, ■ ■ ■ , fn (numite si regresări) • se dau date {măsurători) (si,yi), i = 1, ■ ■ ■ , m Formularea problemei: găsiți coeficienții xi, • ■ ■ ,xn G R: xifi(s,)H xnfn{sj) y, \/i = l, -,m i e găsiți combinații liniare de funcții care fit datele CMMP fit: alegem x care minimizează eroarea de fitting patratica m (X1 H -H Xnfn{Si) - yi)2 i=l in notatie matriceala: Ajj = f/(s7) problema de optimizare min IIAx — y||2 xeR" daca m > n si A are rang n soluția x* = (At4)-1 At y funcția corespunzătoare este: ^mmpfit(s) — xn aplicații: interpolare, netezire de date, identificare de sisteme, • fit un polinom de grad 0 Formularea problemei: estimarea poziției obiectului x si estimarea incertitudinii (zona) • pozițiile ancorelor X1 X2 X3 cunoscute, dar distantele la obiectul necunoscut R1 R2 R3 nu se știu exact: min ll^llp s t : ||x — x/|| reducerea zgomotului echivalenta cu reducerea variației totale Recuperarea perfecta a semnalului original este Concluzie: căutăm “cea mai buna” aproximare a imaginii coru pte Aproximarea se determina utilizând algoritmi de optimizare (pentru niveluri rezonabile de zgomot avem recuperare “aproape” perfecta) Modelare pentru semnale 2D (tip imagine): Asociem imaginii o matrice X G Rmxn sau un vector x G Rmn Dimensiunile (m, n) reprezintă numărul de pixeli pe linii si coloane Ex : Pentru o imagine cu dimensiuni 640 x 480 pixeli asociem X g R640x480 Semnale cu detalii excesive au variație totala mare (i e integrala gradientului absolut al semnalului este mare) reducerea zgomotului echivalenta cu reducerea variației totale: min ||| Y - X||2 + TV(Y) Ye^mxn 2 TV = “total variation" funcție cuantificare variației totala = Ely/+i,j - yij\ + |y/,j+i - x/jl; ij Rezultate (noisy image/denoised image): O imagine se reprezintă numeric printr-o matrice 1 Aproximare matrice cu matrice de rang mai mic via DVS: fie A G Rmxn cu rang(A) = r si A = UY VT = căutăm matricea A avand rangul p n — p Apoi: dim Spân {vi, • • • , vp+i} = p + 1 Cele 2 spatii se intersectează: 3z : ||z|| = 1, Bz = 0, z G Span{vi, ■ ■ ■ , vp+i} p+1 p+1 (A - B)z = Az = ^2 ^iUiV^z = z)uî (z ± vp+2) /=1 i=l p+1 Cum u^uj = 6ij avem 11(A-B)z112 = £ct2(v,rz)211Ui112 > a2p+111z112 i=i M - B\\ = „max ||(A - B)x\\ > ||(A - B)z|| " > ap+1 = ||A - Â|| Interpretare: pentru A = ULVr = 'Ș2'i=1criUivȚ componentele principale de rang 1, UjV^, sunt ordonate in funcție de importanta lor data de at Fie imaginea originala 359x371 pixeli (Melancolia de Durer, 1514): pe care o putem scrie ca o matrice A de dimensiune 359 x 371 (fiecare intrare (/,j) reprezintă intensitatea de gri a pixelului (/,j) cu valori intre 0 (negru) si 255 (alb)), ce poate fi descompusa via DVS: A = LTLVt, unde U este 359 x 359, E este 359 x 371 si V este 371 x 371 Imaginea originala are 359x371 = 133 189 pixeli Dar matricea A poate fi scrisa ca suma de componente principale: A = (TlUivT + ff2U2^2 + H anunv^, unde fiecare matrice utvT (componenta principala) de rang 1 este de marimea matricei originale A Dar pentru ca valorile singulare cr2 > • • • > crn > 0, este posibila o compresie semnificativa a imaginii, atata timp cat spectrul valorilor singulare are doar cateva intrări mari si restul sunt mici! Spectrul lui A conține doar 100 — 200 componente principale max l Putem deci reconstrui fidel imaginea folosind doar o submultime de componente principale De exemplu, putem recupera o imagine asemenatoare cu doar de componente principale folosind comenzile din Matlab: [U,S, V] = svd(A) si B = U(:,l : 20) S(1 : 20,1 : 20) V(:,l : 20) Imaginea B folosește doar din memoria imaginii A\ 20 x 359 + 20 x 371 + 20 = 14 620 vs 359 x 371 = 133 189 pixeli Putem recupera o imagine fidela cu doar de componente principale, in loc de 359 cat are imaginea originala (i e 50% informație din imaginea originala): [t/,5, V] = svd(/\) si B= U(:, 1 : 200)5(1 : 200,1 : 200) V(:,l : 200) Simplificare model: fie modelul y = Ax + 1/ A G R10°x30 are vs a 10> 7> 2, 0 5, 0 01, ■ ■ ■ , 0 0001 ||x|| 5 • modelul simplificat: y = a/U/v^x + v (i e descris de o matrice de rang 4, nu rang 30!) Alta interpretare pentru a,-: min IIA — B\\ B:rang(B) Rp: min ||X||0 s l A(X) = b => min ||X||* s l X(X) = b X^KmXn X£RmXn relaxare convexa • Restricted Isometry Property de ordin r: cea mai mica constanta 5r{A) a i pt orice matrice X de rang cel mult r avem: 1 - 6r(A) 1 Atunci unica soluție a problemei originale neconvexe este matricea de rang cel mult r satisfacand A(X} = b • Presupunem ca r > 1 satisface d$r - 0, si pozitiv semidefinita /IM In cadrul cursului vom utiliza următoarele notatii: Funcții cu litere mici: f, g, h: Rn —> R Gradientul, respectiv matricea Hessiana a funcției continuu diferentiabile f: Rn -G R: Vf(x) G Rn, V2f(x) G Rnxn Hessiana este matrice simetrica (i e in Sn)l dxi 92xi Vf(x) = df(x) 9xn _ , V2f(x) = d2f(x) _dxndxi d2fM~ dx\dxn 92f(x) d2Xn _ Teorema de medie: fie funcția g : R —> R atunci exista 0 G [a, b] a i g(b) - g(a) = g'^b- a) = // g'(r)dT => pt f : R" -> R f (y) = f (x) + (Vf (x + 6»(y - x)), y - x) 0 G f(y) = f(x)+/ (Vf(x + r(y-x)),y-x)dr f* = min f(x) xeR" s L: gi(x) R constrângeri de inegalitate si egalitate g,-, hj: Rn —> R Mulțimea fezabila X = {x : g/(x) 0 • Exemplu de acoperire convexa: pentru o mulțime formata din trei puncte S = {x, y, z} avem Conv(S) = aj = 1 iEl,T finit Hiperplan si semiplan: mulțimi convexe definite astfel hiperplan: {x GR" : aTx = b} , a 0, a G Rn, b G R semiplan: {xGR" : aTx > b] sau {xGR" : aTx b aT x a = r & b = 1 Poliedru: • mulțime convexa definita de p hiperplane sau m semiplane |x G Rn : aȚx = bj \fi = 1, , p, c? x o, fij > o v/,j| • Politop: poliedru mărginit, e g simplex : {x 6 R3 : X1+X2+X3 0} Bila: mulțime convexa definita de o norma || • ||, un centru xc si o raza r: B(xc,r) = {xGR": ||x - xc|| - 0 {xeR": (x — xc)TQ_1(x — xc) 0 avem ax e K Con convex: K este con si mulțime convexa Exemplu con Lorentz- Cn = {[xr t]r E R"+1 : || Exemplu conul Lorentz pentru n = 2: 1 -1 Acoperire conica: Con(S) = {E/eI,z finit :x, GS,a( >0 Conul dual: K* = {y : (x,y) > 0 Vx G K} Daca K = K* atunci K este con auto-dual Exemple conuri: Mulțimea Rn este un con, iar conul sau dual este (Rn)* = {0} R^_ = {x 6 R" : x > 0} se numește conul orthant si este auto-dual in raport cu produsul scalar uzual (x,y) =xTy, i e (R" )* = R" Con pozitiv semidefinit: mulțime formata din matrice pozitiv semidefinite, notatie Sț = {X G Mnxn : X = XT, X 0} X G Sț xTXx >0 VxgR” X,Y G S+,a G xt(qX+(1—n) Y)x = nxTXx+(l—o)xtYx > 0 • Exemplu: y y G S2 o x, z > 0 & xz - y2 > 0 • Observație: conul S” este auto-dual, i e S” = (S”)* in raport cu produsul scalar (X, Y) = Trace(XV) Operații ce conserva convexitatea: intersecția: fiind data o familie de mulțimi convexe atunci P|/e/ S; este o mulțime convexa suma si scalarea: fiind date doua mulțimi Si,S2, atunci Si + S2 = {x + y : x G Si,y G S2} este o mulțime convexa Mai mult, aS = {ax : x G S} este convexa daca S este convexa si a G R fie o funcție afina f(x) = Ax+ b si S mulțime convexa, atunci f(S) = {f(x) : x G S}, este convexa Similar, preimaginea: f_1(S) = {x : f(x) G S} este si ea convexa Exemplu intersecție mulțimi convexe: Inegalitate matriceala liniara (LMI): Consideram o funcție G : Rm —> S/p, G(x) = A) + 12T=-vxi^i> unde x-, sunt componentele unui vector x G Rm, iar matricele A), ■ ■ ■ ,Am G Sn sunt simetrice Expresia G(x) G 0 se numește inegalitate matriceala liniara Mulțimea {x G Rm : G(x)G 0} este o mulțime convexa, fiind o preimagine a lui S™ prin G(x) Exemplu: avem inegalitatea ||AI 0 si variabila A Inegalitatea este echivalenta cu Amax(^rA) b oricare ar fi x e Si si aTx aTXQ oricare ar fi x e S Domeniu efectiv: notam R = RU {+oo} O funcție scalara f : Rn —> R are domeniul efectiv descris de mulțimea: domf = {xGR": f(x) 0 astfel incat f(axi + (1 - a)x2) R, f(x) = x2 f(Yl=V2 Alte exemple de funcții convexe Funcția definita de orice norma este convexa, i e f(x) = ||x|| este convexa pe Rn (folosiți definiția) Funcția f(x) = — log(x) este convexa pe R++ Funcția f(x) = max{xi, ,xn} este convexa pe Rn Funcția f(X) = — logdet(X) este convexa pe spațiul matricelor pozitiv definite S^_+ (dificil de aratatl) • Ineg alitatea lui Jensen este o generalizare a convexitatii: f este convexa daca si numai daca domf este convexa si (p ) p aixj R este continuu diferentiabila si domf este o mulțime convexa Atunci, f este convexa daca si numai daca: f(x2) > f(xi) + Vf(xi)r(x2 — *i) Vxi,x2Gdomf ffe) > /(a-’) ~ - re) Condiții de convexitate pentru funcții continuu diferentiabile Condiții de ordinul II: daca f : Rn —> R este de doua ori continuu diferentiabila si domf este mulțime convexa Atunci f este convexa daca si numai daca pentru orice x G domf matricea Hessiana este pozitiv semidefinita, adica: V2f(x) M Vx G domf Exemplu condiții de ordin I si II: Functia f (x) = — log(x) este convexa pe R++ = {xeR: x > 0} deoarece V2f(x) = ^ > 0 oricare ar fi x > 0 Functia patratica f(x) = ^xTQx + qTx + r este convexa pe Rn daca si numai daca Q 0, deoarece pentru orice x gR" Hessiana V2f(x) = Q (daca Q este simetrica) Se observa ca orice functie afina este convexa si, de asemenea, concava Condiția domf mulțime convexa impusa in toate teoremele precedente este necesara De exemplu, consideram functia f(x) = 1/x2 avand domf = R \ {0} mulțime neconvexa Observam ca Hessiana V2f(x) = 0, dar functia nu este convexa Operații ce conserva convexitatea funcțiilor: Daca fi si sunt funcții convexe si 01,02 > 0 atunci oi fi + 02 f2 este de asemenea convexa Daca f este convexa atunci g(x) = f(Ax + t>) (adica compunerea unei funcții convexe cu o funcție afina) este de asemenea convexa: |||Ax - b||2 = ^xtAtAx - bTAx - ^bTb Fie f : Rn x Rm —> R astfel incât funcția f(-,y) este convexa pentru orice y G S C Rm Atunci următoarea funcție este convexa: g(x) = sup f(x,y) yes Compunerea cu o funcție convexa monotona unidimensionala: daca f : Rn —> R este convexa si g : R —> R este convexa si monoton crescătoare, atunci funcția g o f : Rn —> R este de asemenea convexa • Daca domf convex si exista o constanta a > 0 a i f(axi + (1 - a)x2) f(x) + (Vf(x),y-x) + |||x - y||2 Vx,y Lema 2: In cazul funcțiilor de doua ori diferentiabile (i e f G C2), relația de convexitate tare este echivalenta cu: V2f(x) aln Vx G domf Daca f este tare convexa atunci: 2 Daca domf = Rn, atunci f are un punct de minim global unic: î||x-x-||2 R o funcție patratica, i e f(x) = |xrQx + (q,x) Observam expresia Hessianei In concluzie, pentru funcțiile patratice, cu Q >- 0, constanta de convexitate tare este: q2 > 0, atunci constanta de convexitate tare este a = q? Fie funcția f : Rn —> R, atunci funcția conjugata, notata cu f*, se definește prin f*(y)= max yTx-f(x) xEdom r x ✓ F{x,y) Funcția f* este convexa indiferent de proprietățile lui f Exemple: Pentru funcția patratica convexa f(x) = ^xTQx, unde Q 0, funcția conjugata are expresia: f*(y) = Pentru funcția f(x) = — logx, conjugata sa este data de expresia: r*/ x ( x f-l-log(-y) daca y o I oc altfel Ion Necoara 2014 Optimizarea matematica = = Urcarea unui munte God's Algorithm = Teleportation algoritmi numerici de optimizare: iterativi Pentru problemele prezentate in cursurile I si II vom studia algoritmi numerici (iterativi) de optimizare Un algoritm de optimizare are următoarea structura: Alege un punct intial xq Pasul k\ avand xk, se calculează un vector dk, numit direcție, un scalar ak, numit lungimea pasului, si se actualizeaza iterația: xk+1 = xk + akdk Daca satisface o anumita condiție, atunci ne oprim • Direcția dk se obține diferit, specific fiecărui algoritm • Lungimea pasului insa se alege astfel incat f(xk+i) sa descreasca cat mai mult posibil fata de f(xk), i e ak = arg min f(xk + adk) Astfel, studiem metode de cautare unidimensionala, i e pentru funcții obiectiv-unidimensionale f : R —> R min f(a) agR Pentru punctul de optim folosim notatia a* = arg min f(a) agR Algoritmi de optimizare unidimensionala se bazeaza pe una din următoarele strategii: Cautare directa Interpolare: aproximarea lui f(x) printr-un polinom folosind valorile funcției si/sau derivatele ei in anumite puncte Cautarea directa: este bazata pe următorii pași Se identifica un interval inițial [ao, M care sa il conțină pe a* Intervalul inițial se reduce iterativ [a/j, bk] —> [a/j+i, b/ 0 daca f(ao + ho) f(ao), atunci ne deplasam din oo inapoi pana când valoarea funcției creste Metodele de cautare unidimensionale se bazeaza pe unimodalitatea funcției: Daca exista a* G [a, b] astfel incat f este strict descrescătoare pe intervalul [a, «*] si strict crescătoare pe intervalul [a*, b], atunci f se numește funcție unimodala pe intervalul [a, b] Teorema: fie f unimodala pe [a, b] si fie doi scalari 01,02 G [a, b], cu Oi ffe), atunci este interval de unimodalitate pentru f Folosing unimodalitatea, metoda forward-backward are următorii pași (presupunem ca a* G R+) Pas 1 fie un oq G [0, oo), ho > 0 si coeficientul multiplicativ t > 1 (adesea se alege t = 2) Evaluam f(ao) = fo si k = 0; Pas 2 Comparam valorile funcției obiectiv Actualizam ak+1 = ak + hk si evaluam fk+1 = f(ak+1) Daca fk+1 atunci stop si returneaza [a, b]; altfel a = r], rj = b, sărim la Pas 2; Pas 4 Pas backward: rj = otk, b = ak+i, hk+i = thk, a = b — hk+i Daca f(a) > atunci stop si returnam [a, b]; altfel, b = r], si rj = a, sărim la Pas 2 Exemplu metoda forward-backward pentru f (a) = s/n(27ra) + cos(57ra) 0 2 0 4 0 6 0 8 1 1 2 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha alpha alphak=0 400000 alphakt1 =0 400100 k=0 h=0 000100 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha alpha alphak=0 400100 al phakt1 =0 400300 k=1 h=0 000200 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha alpha alphak=0 400300 al phakt1 =0 400700 k=2 h=0000400 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha alpha alphak=0 400700 alphakt1 =0 401500 k=3 h=0 000800 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha alpha alphak=0 401500 alphakt1 =0 403100 k=4 h=0001600 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha alpha alphak=0 403100 al phakt1 =0 406300 k=5 h=0 003200 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha alpha alphak=0 406300 alphakt1 =0 412700 k=6 h=0006400 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha alpha alphak=0 412700 al phakt1 =0 425500 k=7 h=0 012800 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha alpha alphak=0 425500 alphakt1 =0 451100 k=8 h=0025600 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha alphak=0 451100 al phakt1 =0 502300 k=9 h=0051200 alpha 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alpha al phak=0 502300 al phakt1 =0 604700 k=10 h=0102400 I f(x)| alpha 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) alphak=0 604700 al phakt1 =0 809500 k=11 h=0204800 alpha alpha 0 2 0 4 0 6 0 8 1 1 2 1 Exemplu metoda forward-backward pentru f(a) = s/n(27ro!) + cos(57ra) a f i] La iterația k, metoda secțiunii de aur determina intervalul [a/(_i_i, b/ f(fJ-k) atunci ak+i = Xk si bk+i = bk- întrebare: cum alegem acum A/j si [ik ? Răspuns: impunem condițiile distantele de la Xk si respectiv /ik la capetele intervalului [a/t, bk] sunt egale: bk ~ = b-k ~ 3k- (1) rata de micșorare a lungimii intervalelor de incertitudine la fiecare iterație este aceeași, rezultând bk+i - ak+i = r(bk - ak), unde r G (0, 1) (2) este necesara o singura evaluare a funcției obiectiv pentru o noua iterație Condiția 1: distantele de la Xk si respectiv p,k la capetele intervalului [ak, bk] sunt egale: bk - Xk = P-k - ak■ Reamintim cazul (i) din unimodalitate: f(Xk) 0 rezulta t = v/52 1 — 0-618 (numărul de aur) Astfel, metoda secțiunii de aur consta in: Pas 1 Determinam [a1, b1], alegem precizie 6 Calculam X1 = a1 + 0-382(b1 — a1), b1 = a1 + 0-618(b1 — a1), evaluam f(X1) si f(b1), initializam k = 1 Pas 2 Daca f (Xk) > f (bk) Pas 3; altfel Pas 4 Pas 3 Daca bk — Xk fkj , iar din iterațiile l-k-J+l metodei secțiunii de aur se obține: bj+i - aj+1 = Fk J (bj - aj) ’k-j+l Impunem o toleranta asupra lungimii Luând in considerare: bk~3k — -p-(bk-i-ak-i) — -p p- ■ ■ t~2 ~3 avem Fk > bl~ai => se pot calcula d Asimptotic insa lim^oo = t (= (valabil in ambele cazuri) intervalului final: bk — ak ak+i = ak - ——-r"(a/ ak+i = otk~ (a/J f (a/c-i) — ' («/ 0 Exemplu de pas metoda falsei poziții: Xk+1 Xk-1 x Fiecare student alege o aplicație a i sa se modeleze ca o probleme de optimizare unidimensionala Cerințe: descrieți in detaliu aplicația aleasa dati formularea matematica ca o problema de optimizare utilizați cel puțin 2 algorithmi de optimizare unidimensionala dati la curs pentru rezolvarea ei dati codul Matlab pentru cei 2 algoritmi comparați rezultatele cu funcția Matlab fminbnd si fminunc comentați rezultatele obținute (grafice, tabele) Ion Necoara 2014 Consideram problema de optimizare: (UNLP) : min f(x) xeR" Presupuneri: f eC2 domf mulțime deschisa in Rn Problema de interes consta in gasirea punctelor de minim din domf Reamintim ca un punct x* se numește punct de minim: global f(x*) 0 a i f(x*) 0 suficient de mic a i din continuitatea gradientului avem: Vf(x + rd)Td contrazice ipoteza inițiala Datorita continuității gradientului, putem alege un t > 0 suficient de mic a i Vr G avem: Vf(x* + rd)Td = -Vf(x* - rVf(x*))rVf(x*) 0 si X2 = — 2xi obținem 7(xi,X2) 0 = 7(0,0) Concluzie: (0,0) punct de inflexiune Ce tip sunt restul punctelor staționare? Cum putem analiza natura punctelor in cazul problemelor de mari dimensiuni? După cum am văzut anterior, natura punctelor staționare nu poate fi stabilita utilizând doar informație de ordinul I pentru a stabili ca un punct staționar este punct de extrem (minim sau maxim) este necesara informație despre Hessiana Teorema (Condiții necesare de ordinul II): Fie f G C2 si x* G domf un punct de minim local al problemei (UNLP) Atunci avem: V2f(x*) > 0 Demonstrație: Presupunem prin absurd ca pentru un x* punct de minim local avem ca V2f(x*) > 0 nu este satisfăcuta si drept urmare 3d a i drV2f(x*)d 3t > 0 a i Vr G avem: drV2f(x* +rcT)d 30 G a i : f(x* + td) = f(X^3-t\7f(x*)Td +^t2 dT\72f(x* + 0d)d - 0, atunci x* este punct strict de minim local al problemei (UNLP) Demonstrație: Notând cu Amin > 0 valoarea proprie minima a Hessianei V2f(x*) avem: drV2f(x*)d > Amin||d||2 VdGR" Folosind in continuare aproximarea Taylor putem scrie: f(x* + d) - f(x*) = Vf(x*)rd + |drV2f(x*)d + TI (||d||2) > + «(!! o pentru orice d G Rn suficient de mic (reamintim ca avem relația lim = 0) Deci x* este punct de minim local Wl 7 Ca o consecința a condițiilor suficiente de ordinul II, putem stabili natura unui punct staționar x* in functie de proprietățile Hessianei, astfel: punct de minim local V2f(x*) >- 0 punct de maxim local V2f(x*) - f(x*) + Vf(x*)r(x - x*) = f(x*) Vx G Rn =o Concluzie: orice punct de minim local al unei probleme convexe este punct de minim global! Probleme convexe sunt mult mai ușor de rezolvat! Consideram problema (UNLP) convexa: min f(x) xeR2 — X1 + — Xl*2 + 2X2 In continuare obținem ca Vf(x) si V2f(x) au următoarele expresii: Vf(x) = 4xx + 4xi - x2 -xi + x2 si V2f(x) = 12x2 + 4 -1 -1 1 Observam ca V2f(x) >- 0 VxgR2 deci f este convexa Conform condițiilor suficiente de ordinul I, soluțiile sistemului Vf(x) = 0 sunt puncte de minim global ale problemei (UNLP) (0,0) punct de minim global Ion Necoara 2014 Problema de optimizare neliniara fara constrângeri (UNLP) : f* = min f(x) xeR" Condițiile de ordinul I: Vf(x) = 0 (n ecuații cu n necunoscute) Exemplu (QP): min 0 5xrQx — qTx, xeR" Q matrice inversabila Vf(x) = Qx — q = 0 soluție unica x* = Q-'q valoarea optima f* = — 0 5qTQ_1q (UNLP) : f* = min f(x) xeR" Sistemul de n ecuații cu n necunoscute Vf(x) = 0 nu poate fi rezolvat analitic! metode iterative de rezolvare Algoritm iterativ: generează un sir (x/()/(>o; fiecare punct calculat pe baza punctelor anterioare: xk+1 = M(x0, ■■■ ,xk) algoritm convergent global: pentru orice x0 șirul xk converge la o soluție x* algoritm convergent local: pentru orice xq apropiat de o soluție x* șirul xk converge la x* (UNLP) : f* = min f(x) xeR" Algoritm iterativ: generează un sir de puncte (x/ o; fiecare punct calculat pe baza punctelor anterioare: xk+1 = M(x() ■■■ ,xk) Algoritmul numeric colectează informație printr-un oracol si manipulează răspunsurile oracolului oracol de ordin zero O®, furnizează informație bazata pe evaluarea funcției obiectiv f(x) oracol de ordin intai C\: furnizează informație bazata pe evaluarea funcției obiectiv f(x) si gradientului Vf(x) oracol de ordin doi C>2'- furnizează informație bazata pe evaluarea funcției obiectiv f(x), gradient Vf(x) si Hessiana V2f(x) Algoritm de optimizare iterativ gaseste o soluție aproximativa cu o acuratețe prestabilita e acuratețea e reprezintă si criteriul de oprire al algoritmului criteriul de oprire 1: ||Vf(x/f)|| 0 ales la pasul inițial k = 0 • Complexitatea analitica a algoritmului de optimizare iterativ: numărul total de apeluri ale oracolului pentru gasirea unei e-solutii • Rata de convergenta: viteza cu care șirul xk converge la soluția x*; ordinul de convergenta q definit ca: Hxk+i - x*|| & (3\\xk - x*\\q convergenta subliniară: \\xk — x*|| 0 ||x/(_i_i — x*|| & (3\\xk — x*\\q • distanta se reduce cu q zecimale convergenta subliniară: *k = j; => (1; 0-5; 0 33; 0 25; 0 2) conv liniara: x/ (0 33; 0 11; 0 03; 0 01; 0 004) conv superliniara: Xk = => (1; 0 16; 0 05; 0 01; 0 002) conv patratica: xk=^k =>(0 1; 0 01; 0 001; 0 0002; 0 00001) Definiție Se da: spațiul metric (X p), submultimea S C X si metoda JVt de tip aplicație punct-multime (JVt : X —> 2X) Definim funcția descrescătoare : X —> R pentru perechea (S, JVt): Vx G S si y G JVt(x) avem (y) (x) Vx S si y G JVt(x) avem (y) (x) Exemplu: pentru problema minxf(x), definim S = {x : Vf(x) = 0} (mulțimea punctelor staționare) si = f (de obicei o metoda de optimizare alege x/j+i a i f(x/(+i) 2X este inchisa in xo daca pentru orice Xk —> xo si yk —> yo cu yk G M(xk) avem yo G JVt(xo) M este inchisa daca este inchisa in toate punctele din X Exemplu: aplicația punct-multime x/j+1 G [—\xk\, |x^|] este inchisa Theorema de convergenta generala Fir șirul x/j+i G M(xk) iar S mulțimea soluțiilor, satisfacand: Xk se afla intr-o mulțime compacta M este aplicație punct-multime inchisa exista o functie continua descrescătoare pentru Atunci punctele limita ale șirului Xk aparțin mulțimii soluțiilor S xk+1 = xk + a kdk, unde dk este direcție de descreștere pentru f in xk, i e Vf(xk)Tdk 0 f(xk + otdk) metoda Wolfe: (Wl) : f(xk +akdk) c2Vf(xkyrdk cu ci O si p, ci G (O, 1) cat timp f(xk + adk) > f{xk) + cia/fVf(x/f)rd/f a = pa Fie o funcție continuu differentiabila f (i e f G C1), atunci gradientul Vf este continuu Lipschitz cu parametrul L > 0 daca: ||Vf(x) - Vf(y)|| ||V2f(x)|| ^-||x-y||2 + ^-||Vf(y)-Vf(x)||2 Fie f : Rn —> R o functie patratica, i e f(x) = |xrQx + (q,x) Observam expresia gradientului Aproximam constanta Lipschitz a funcției f: \\Qx + q — Qy — q\\ = ||Q(x - y)|| Q2 > 0, atunci constanta Lipscithz L = R definita de f(x) = log (1 + eaTx^ Observam expresia gradientului si a matricii Hessiene Pentru orice constanta pozitiva c > 0 avem ^^2 (c2 - l)Vf(xk)Tdk Aplicând inegalitatea Cauchy-Schwartz obținem: ||Vf(xk+i) - Vf(xk)|| ■ \\dk\\ > (c2 - l)Vf(xk)Tdk Din proprietatea de Lipschitz gradient avem: ||Vf(x/(+i) — Vf(x/()|| (c2 - l)^f(xk)Tdk i e c2 - l^f{xk)Tdk a‘- L ||4||2 ■ Din condiția Wolfe (Wl) avem: fi \ oo Ion Necoara 2014 Informația ce indica comportamentul unei funcții f gR" ->R intr-un punct xGR" se poate clasifica: Informație de ordin 0: f(x) Informație de ordin 1: f(x), Vf(x) Informație de ordin 2: f(x), Vf(x), V2f(x) Fie algoritmul iterativ definit de x/j+i = in functie de ordinul informației utilizate in expresia lui jVt: Metode de ordin 0: f(x/() Metode de ordin 1: Vf(x/ O(^p-) (sublinear - gradient Lipschitz) • ^((i^)k) -> O((l — a/t)^) (liniar ’ tare convex + grad- Lip ) Fie funcția f : Rn —> R diferentiabila (UNLP) : min f(x) xeR" Iterație Metoda Gradient: Complexitate pe iterație O(n) daca evaluarea Vf(x) este ieftina! Metoda gradient rezolva probleme de IO9 variabile din motoare de cautare, procesarea de imagine ( Interpretare: metoda de descreștere cu direcția d = — Vf(xk), deci f(x/t+i) 0 daca: ||Vf(x) - Vf(y)|| ^-||x-y||2 + ^-||Vf(y)-Vf(x)||2 Fie f : Rn —> R o functie patratica, i e f(x) = |xrQx + (q,x) Observam expresia gradientului Aproximam constanta Lipschitz a funcției f: \\Qx + q — Qy — q\\ = ||Q(x - y)|| R definita de f(x) = log (1 + eaTx^ Observam expresia gradientului si a matricii Hessiene Pentru orice constanta pozitiva c > 0 avem ^^2 0) si mărginită inferior Daca alegem lungimea pasului ak astfel incat satisface condițiile Wolfe, atunci șirul Xk generat de metoda gradient satisface: lim Vf(x/() = 0 / oc Demonstrație: Se observa ca unghiul gradientului fata de direcția metodei gradient (antigradient) este dat de 0k = n Din Teorema de Convergenta Globala a metodele de descreștere avem: ^cos2^||Vf(x/ 0 k>0 Rezulta: Vf(x/() —> 0 când k —> oo Teorema 5: Fie f diferentiabila cu Vf Lipschitz continuu (constanta Lipschitz L > 0) Alegem lungimea pasului = ț Atunci rata de convergenta globala a șirului Xk generat de metoda gradient este subliniară, data de: min ||Vf(x/)|| 0) Exista un punct de minim local x*, astfel incat Hessiana in acest punct satisface a In = 0 Teorema 7 Fie f functie , diferentiabila cu Vf Lipschitz continuu (constanta Lipschitz L > 0) Daca alegem lungimea pasului constanta ak = ț, atunci rata de convergenta globala a șirului Xk generat de metoda gradient este subliniară, data de: f(xk) -f* 0, atunci rata de convergenta globala a șirului Xk generat de metoda gradient cu pas | este liniara, data de: \\xk -x*||2 ||x-y||2 + ^—||Vf(x)-Vf(y)||2 a + L a + L Aceasta relație conduce la: llx^+1 — X*||2 = \\xk - l/Z-V/Xx/J — X* ||2 = \\xk - x* ||2 - 2/L(Vf(xk),xk-x*) + 1/L21| Vf(xk)||2 O’Ă) l|x‘-x*l|2+(Ă-Z(Ăz))l|vfWI12 coerciva+Vf(x* )=0 s -v z Pe de alta parte, in valoarea funcției (cu gradient Lischitz) avem: - vil2 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 ( = |(0-5*i2 +7*2)) , cu 7 > 1 min f(x) xeR2 Funcția f tare convexa (a = 0 5) si gradient Lipsctitz (/ = 7) convergenta liniara Metoda Gradient cu pas constant a = ț Punct inițial xq = [|7 1], 7 = 5/3 min f(x) I = -(0 5x? + rx|) ) xeR2 \ 2 2 J Metoda Gradient cu pas ideal ak = argmin f(xk — a\7f(xk)) a>0 Punct inițial xq = [|7 1], 7 = 5/3 min f(x) I = -(0 5x? + rx|) ) xeR2 \ 2 2 J Metoda Gradient cu pas ideal ak = argmin f(xk — a\7f(xk)) a>0 Punct inițial xq = [|7 1], 7 = 5/3 min f(x) I = -(0 5x? + rx|) ) xeR2 \ 2 2 J Metoda Gradient cu pas ideal ak = argmin f(xk — a\7f(xk)) a>0 Punct inițial xq = [|7 1], 7 = 5/3 min f(x) I = -(0 5x? + rx|) ) xeR2 \ 2 2 J Metoda Gradient cu pas ideal ak = argmin f(xk — a\7f(xk)) a>0 Punct inițial xq = , 7 = 5/3 min f(x) I = -(0 5x? + rx|) ) xeR2 \ 2 2 J Metoda Gradient cu pas ideal ak = argmin f(xk — a\7f(xk)) a>0 Punct inițial xq = [|7 1], 7 = 5/3 min f(x) I = -(0 5x? + rx|) ) xeR2 \ 2 2 J Metoda Gradient cu pas ideal ak = argmin f(xk — a\7f(xk)) a>0 Punct inițial xq = [|7 1], 7 = 5/3 xGR2 exi+3x2-0 1 _|_ exi-3x2-0 1 + e-*1-0'1 funcție obiectiv convexa (nu este tare convexa, nu are gradient Lipschitz pe R2) Metoda Gradient cu pas backtraking / ideal backtracking line search exact line search Rata de convergenta slaba a metodei gradient reprezintă motivația dezvoltării de metode cu performante superioare Metoda de Gradient Accelerat (Nesterov 1983) - cu un ordin mai rapida decât gradientul clasic in cazul problemelor convexe Metoda de Gradienti Conjugați (Lanczos, Hestenes, Stiefel 1952) - pentru QP convex soluția in n iterații Metoda de Gradient Accelerat *k+i = Yk~ f(yk) yk+i = xk+1 + /3k(xk+1 - xk) unde se iau punctele inițiale xq = yo si fik ales in mod adecvat: e g sub convexitate tare putem alege fik = yi-+>/â Observație: Costul iterației similar cu cel al metodei gradient clasice! Teorema 8 Fie f o funcție , diferentiabila cu Vf Lipschitz continuu (constanta Lipschitz L > 0) Rata de convergenta globala a șirului xk generat de metoda gradient accelerat este subliniară, data de: f(*k) -f* 0, atunci rata de convergenta este liniara, data de: k f(xk) -f* (sublinear - gradient Lipschitz) • > 0((1 — a/t)^) G'n'ar ’ tare convex + grad Lip ) O(^) O(^) versus O((^)/ 1 Funcția obiectiv f tare convexa (a = 0 5) si gradient Lipschitz (/ = 7) Metoda Gradient cu pas constant = 1/L Metoda Gradient Accelerat cu pas constant ak = 1/L si fik = (VL - y/a)/(yT + ^/ă) Ambele metode de ordinul I converg liniar Punct inițial xq = [|7 1], 7 = 5/3 Metoda Gradient cu pas constant ak = 1/L Metoda Gradient Accelerat cu pas constant ak = 1/L si Pk = (vT - v/ă)/(vT + v/ă) Metoda Gradient cu pas constant ak = 1/L Metoda Gradient Accelerat cu pas constant ak = 1/L si Pk = (vT - v/ă)/(vT + v/ă) Metoda Gradient cu pas constant ak = 1/L Metoda Gradient Accelerat cu pas constant ak = 1/L si Pk = (vT - v/ă)/(vT + v/ă) Metoda Gradient cu pas constant ak = 1/L Metoda Gradient Accelerat cu pas constant ak = 1/L si Pk = (vT - v/ă)/(vT + v/ă) Metoda Gradient cu pas constant ak = 1/L Metoda Gradient Accelerat cu pas constant ak = 1/L si Pk = (vT - v/ă)/(vT + v/ă) Metoda Gradient cu pas constant ak = 1/L Metoda Gradient Accelerat cu pas constant ak = 1/L si Pk = (vT - v/ă)/(vT + v/ă) Metoda Gradient cu pas constant ak = 1/L Metoda Gradient Accelerat cu pas constant ak = 1/L si Pk = (vT - v/ă)/(vT + v/ă) Metoda Gradient cu pas constant ak = 1/L Metoda Gradient Accelerat cu pas constant ak = 1/L si Pk = (vT - v/ă)/(vT + v/ă) min f(x) x£Rn m = log Ș2 eaîx+bi j=i Funcția obiectiv are gradient Lipschitz, dar nu e tare convexa Date generate aleator cu m = IO4 si n = IO3 Metoda Gradient si Gradient Accelerat cu pas ak = 1/L Metoda Gradient accelerat nu este o metoda de descreștere! Fie problema de optimizare patratica neconstransa min ^-xtQx — (q,x) xeR"2 x 7 Q este matrice simetrica, pozitiv definita: soluția optima x* = Q_1q • rezolvarea sistemului Qx = q Definim residul r = Qx — q Doi vectori di si c/2 sunt daca Qdz = 0 Daca {di, - ,dn} sunt Q-ortogonali (deci liniar independenți), atunci exists scalarii ak a i n X* = '^akdk cu k=l djQx* ak } = Span{/-0, Qro, • • -,Qkr0} d^ Qdi = 0 pentru orice / rk+1 rk rk Algoritmul MG (Se da punctul de start x0 si acuratețea e Se calculează o e-solutie optima pentru problema de optimizare minx f(x) (= 10x® + 30x® + x2 + cu MG-ideal ) 0 function [■] = MG-ideal(xO, e) 1 obj = @(x) 10 * x(l)6 + 30 * x(2)6 + x(l)2 + 50 * x(2)2 2 grad = @(x) 3 x = xO, tg = xO 4 while(norm(grad(x)) > e) 1 objn = @(a) obj(x — agrad(x)) 2 a* = fminbnd(objQ,, 0,1) 3 x = x - a* * grad(x); tg = [tg; x] 5 end while 6 x = —0 2 : 0 1 : 0 2; y = —0 2 : 0 1 : 0 2; [X, V] = meshgrid(x, y) 7 Z = 10 * X6 + 30 * Y6 + X2 + 50 * Y2 8 figure; plot(tg(l,:), tg(2,:)); hold on; contour(X, Y, Z) Ion Necoara 2014 Scurt istoric: antichitate Metoda babiloniana sau metoda lui Heron: dandu-se un număr c E R, aceasta metoda consta in calcularea iterativa a lui x = y/c O x2 — c = 0: 1 xk+i = 2 c A xk H -I xk / = - 2^ R si problema de optimizare aferenta: min f(x) xeR" Din condițiile de optimalitate de ordinul I pentru probleme neconstranse amintim: x* este punct de minim local daca Vf(x*) = 0 si V2f(x*) >- 0, unde V2f(x*) denota Hessiana lui f(x) in punctul x* Astfel, pentru a afla x* trebuie sa mai intai aflam soluția ecuației: F(x) = Vf(x) = 0 • Interpretare: Aplicând metoda lui Newton-Raphson in acest caz, rezulta iterația: Fiind metoda iterativa, metoda lui Newton se poate scrie drept: *k+i = *k + dk unde • Interpretare: daca Vf(x)2 >- 0, atunci (d/j, Vf(%/ - O, atunci q(y) este o funcție patratica strict convexa Din condițiile de optimalitate pentru probleme QP strict convexe, rezulta: Vq(y) = 0 => V f(xk) + V2 f(xk) (y -xk) = 0 Din moment ce y = xk+1, prin rezolvarea ecuației anterioare in y rezulta din nou iterația: xk+1 =xk- (V2f(xk))_1Vf(xk) Observație, din moment ce direcția Newton dk = xk+i — xk, atunci se poate exprima drept: dk = argmin f(xk) + Vf(xk)rd + ^drV2f(xk)d d 2 Din nou, utilizând condițiile de optimalitate, rezulta: c/k = -(V2f(xk))_1Vf(xk) Iterație metoda Newton prin aproximarea patratica Taylor: Observație: metoda Newton pentru funcții f : R —> R este chiar metoda Newton de interpolare din cursul III: xk+1 = argmin q(x) = argmin f(xk) + f'(xk)(x-xk) + -f"(xk)(x-xk)2 = xk-(f\xk)y1f'(xk) Teorema (rata de convergenta locala a metodei Newton): fie f G C2 si x* un minim local ce satisface Vf(x*) = 0 si V2f(x*) >- 0 fie o constanta m > 0 astfel incat: V2f(x*) mln presupunem ca V2f(x) este Lipschitz cu constanta M > 0, i e ||V2f(x) — V2f(y)|| 0 Prima demonstrație data de Kantorovich in 1948! Demonstrație: x* este un minim local —> Vf(x*) = 0 Din teorema lui Taylor in forma integrala avem: Vf(xk) = Vf(x*)+ f Jo \72f(x* + r(xk - x*))(xk - x*)dr Se obține: Xk+1 _x*=xk_x*_ (V2f(xk))-1Vf(x/() = (V2f(xk))-1[V2f(xk)(xk - x*) - Vf(xk) + Vf(x*)] = (V2/r(xk))“1[V2f(x/()(x/(-x*)- f V2/7(x*+t(x/(-x*))(x/(-x*)c/t] Jo = (V2f(xk))_1 f [V2/r(xk)(xk—x*)—V2f(x* + r(x/f-x*))(x/f-x*)cfr] Jo = (V2f(xk))_1 [ [V2f(xk) - V2f(x*+T(xk-x*))](xk-x*)c/T Deoarece ||V2f(jx7c) — V2f(x*)|| = V2f(x*) - M\\xk — x*\\ln h mln - M\\xk — x*\\ln >- 0, sub ipoteza ca \\xk — x*|| 0 Observam ca -y - 0 Din condițiile de optimalitate: Vf(x*) = Qx* + q = 0 => x* = -Q~xq Pornind dintr-un punct xo, si aplicând metoda Newton, rezulta: xi = x0-(V2f(x0)) 1 W(x°) = x0-Q_1(Qx0+q) = -Q~1q = x* deci pentru probleme patratice strict convexe metoda converge intr-un singur pas Dar metoda Newton necesita calculul inversei V2f(x)_1, operație costisitoare O(n3) pentru probleme de dimensiuni mari Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Newton clasic 1 Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Newton clasic Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Exemplu metoda Newton cu pas variabil ales in trei variante, aplicata pe f: R2 —> R, f(x) = exi+3x2_o-1 + exi_3x2_o'1 _|_ e-xi-o i Newton clasic Newton clasic Newton clasic Newton clasic Newton clasic Newton clasic Newton clasic Newton cu backtracking Newton cu backtracking Newton cu backtracking Newton cu backtracking Newton cu backtracking Newton cu backtracking Newton cu backtracking Newton cu backtracking Newton cu backtracking Newton cu backtracking Newton cu backtracking Newton cu pas ideal Newton cu pas ideal Newton cu pas ideal Newton cu pas ideal Newton cu pas ideal 50 40 30 T î 20 10 160 140 120 100 60 40 20 Newton clasic Newton ideal Newton backtracking 0 15 k 20 Newton clasic Newton ideal Newton backtracking 10 15 20 Daca pornim dintr-un punct xq ce nu se afla in vecinătatea lui x*, metoda Newton trebuie modificata pentru a asigura convergenta: Pasul ak se poate alege ideal (i e ak = mina>o f(xk +akdk)) sau prin metoda backtracking Teorema (convergenta globala a metodei Newton): Fie funcția obiectiv f G C2 mărginită inferior si cu gradientul Vf Lipschitz Consideram metoda Newton cu pas variabil ak ales prin backtracking: Xfc+1 = xk- ak(y2f(xk^ 1^f(xk) Presupunem ca Hessiana satisface condiția fi-Jn = complexitate pe iterație O(n) plus costul evaluării Vf(x)! • Metoda Newton cu pas ak se obține prin aproximarea Taylor de ordin II: xk+l =Xk- ak(y2f(xk)Y1'Vf(xk') = argmin f{xk)+\7f{xk)T{y-xk) + ^{y-xk)T\72f{xk){y-xk) y 2ak => complexitate pe iterație O(n3) plus costul evaluării Vf(x) si V2f(x)l min f(x) (= lOxf + 30xo + *1 + 50x|) xeR2 function [] = gradient-Newton-ideal(xO,eps) obj = @(x)10x® + 30*2 + *1 + 50xj; gradient = @(x) ; hessiana = @(x) : 7 7 Metoda Gradient cu pas ideal x = xO; trajectoryg = [x0]; while (norm(gradient(x)) > eps) grad = gradient(x); obja = @(a) obj (x - a grad); a* = fminbnd(obja, 0, 1); x = x - a* gradient(x); trajectoryg = [trajectorygx]; end 7°/, Metoda Newton cu pas ideal x=xO; trajectoryn = [xO] ; while (norm(gradient(x)) > eps) grad = gradient(x); hess = hessiana(x); d-newton = inv(hess)* grad; obja = @(a) obj (x - a *d-newton); a* = fminbnd(obja, 0, 1); x = x - a* d-newton; trajectoryn = [ trajectoryn x] ; end x = -0 2 : 0 1 : 0 2; y = -0 2 : 0 1 : 0 2; [X,Y] = meshgrid(x,y); Z = 10*X "'6 + 30*Y "6 + X "2 + 50*Y ’'2; f igure plot (trajectoryg (1, :), trajectoryg (2,:) , ’ r+- ’, ’ LineWidth ’, 3); hold on plot (trajectoryn (1, :) , trajectoryn (2, : ) , ’ k*— ’, ’ LineWidth ’, 3); legend('Gradient Method’,'Newton Method’); hold on contour(X,Y,Z,’ShowText’,'on','LineWidth',2); end Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Gradient cu pas ideal Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Gradient cu pas ideal Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Gradient cu pas ideal Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| 0 2 0 15 0 1 0 05 o -0 05 -0 1 Gradient cu pas ideal -0 15 -0 2 Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| Comparație intre metoda Newton si metoda gradient, ambele cu pas ideal, pentru f: R2 —> R, f(x) = lOxf + 30xf +x2+50x| 3 ■■ 2- 1- 3 ■■ 2- 1- 3 ■■ 2- 1- 3 ■■ 2- 1- Gradient cu pas ideal 3 ■■ 2- 1- 3 ■■ 2- 1- 3 ■■ 2- 1- 3 ■■ 2- 1- 3 ■■ 2- 1- 3 ■■ 2- 1- 3 ■■ 2- 1- 3 ■■ 2- 1- Gradient cu pas ideal 3 ■■ 2- 1- 3 ■■ 2- 1- Avantajul metodei Newton: convergenta extrem de rapida Dezavantajul metodei Newton: folosește informație de Hessiana V2f(x) si rezolvarea unui sistem liniar (inversa (V2f(x/ 0 atunci Hi >- 0 (i e updatarile BFGS pastreaza pozitiv definitatea) cost pe iterație: O(n2) plus costul calculării Vf(x) BFGS cu pas a = 1 pe cazul scalar coincide cu metoda falsei poziții Teorema (rata de convergenta a metodelor cvasi-Newton): fie x* un punct ce satisface condițiile suficiente de optimalitate de ordinul II presupunem iterația forma x/j+i = xk — Hk este inversabila si satisface condiția Lipschitz: l|Hk(V2f(xk)-V2f(y))|| 0 sau rata liniara daca 7^ > 7 > 0 Consideram problema (UNLP): m min cTx — > log(5/ — aj x) j=i Generam datele aleator pentru dimensiunea n = 100 si m = 500 Newton B=GS 10 10 10 io1 10 1010 10 10 10 io' 10’ 10 io10 10 0 20 40 60 80 100 120 140 cost pe iterație pentru metoda Newton: O(n3) plus costul calculării Vf(x) si V2f(x) cost pe iterație pentru metoda BFGS: O(n2) plus costul calculării Vf(x) metoda gradient: converge cel mult liniar si costul pe iterație este O(n) (operații cu vectori) plus evaluarea gradientului Vf(x) eficient in rezolvarea problemelor de foarte mari dimensiuni n = IO9 metoda Newton: converge patratic si costul pe iterație este O(n3) (rezolvarea unui sistem liniar) plus evaluarea gradientului Vf(x) si Hessianei V2f(x) eficient in rezolvarea problemelor de mici dimensiuni n f diferenția bila xeR" Exemplu: sparse data fitting (aproximarea unui set de date cu funcție liniara): set de date (a,-, b,) a i b,- depinde linear de a,-: b/ = ajx + et \/i = 1 : m cu e,- o eroare datorata procesului de sampling => e g CM MP rar min fT(x) (= ||Ax — b||2 + r||x||i) => A G Rmxn xeR" Aplicații: Magnetic Resonance Imaging (MRI): instrument medical imagistic folosit pt scana anatomia si fiziologia unui organism Image inpainting: technica de reconstrucție imagini degradate Image deblurring: instrument de procesare imagine pt indeparta neclaritatea cauzate de fenomene naturale (mișcare) Genome-Wide Association study (GWA): comparare ADN intre 2 grupuri de persoane (cu/fara o boala), in scopul de a investiga factorii de care depinde o boala min fT(x) (= f(x) + t||x||i) => f diferenția bila xeR" Aproximam f cu o funcție patratica si păstrăm nealterat r||x||i: xk+1=arg rn\nqT(x-xk)(=f(xk)±(^f(xk) x-xk}±h\x-xk\^k±T *lli) Hk = Lln unde L constanta Lipschitz a lui Vf => metoda gradient (accelerat): Nesterov’07 Hk = diag(/_i, ■ ■ ■ , Ln), unde L,- constanta Lipschitz a lui V,-f => metoda gradient pe coordonate: Nesterov’10, Necoara’ll Hk = ^2f(xk) si inlocuieste funcția nediferentiabila ||x||i cu aproximarea diferentiabila (Huber) |x,| o ^/x2 + M2 — /' => metoda Newton: Gondzio’15 Probleme de data fitting cer analiza Big Data (terabytes of data): trilioane de varibile (n) => Cray XC30 MPP supercomputer (ARCHER in Edinburgh este al 25-lea in top 500): 118 080 cores cu performanta 1 642 TFIops/s on LINPACK benchmark ARCHER in Edinburgh: 25-lea in top 500 cu 118 080 cores Comparație metode: gradient pe coordonate paralel (PCDM), gradient accelerat (FISTA), Newton (pdNCG) cu gradient conjugat pt rezolvare aproximativa sistem liniar si quasi-Newton (PSSgb) Ion Necoara 2014 Problemele de estimare si fitting => probleme cu structura, de tipul celor mai mici patrate (CMMP): jVt(x) = r] forma algebrica min -||î?-7Vf(x)||2 x€Rn 2 forma optimizare tj G Rm m măsurători; jVt : Rn —> Rm => modelul considerat; x EE Fi" => parametrii modelului Problema forward: pentru intrări date ale modelului jVt(x) se determina ieșirile corespunzătoare Problema inversa: pentru un vector de ieșiri rj se cauta intrările corespunzătoare, folosind un model ce depinde de setul de parametrii x Aplicații inginerești: aproximare de funcții, identificarea proceselor dinamice, prognoza meteo Problema CMMP isi are rădăcinile in domeniile astronomie si geodezie, apariția sa fiind motivata de aspecte legate de navigarea pe oceane, elementul cheie in dezvoltarea sa reprezentandu-l problema descrierii traiectoriei astrelor pe orbita In 1809 => deter- minarea orbitei planetei “strumf” Ceres => susține ca dispunea de aceasta metoda inca din 1795 Prima expunere concisa a metodei ii aparține lui (1805) 1810 => => teorema de limita centrala ce oferă o corelație intre problema CMMP si distribuția normala Consideram un set de m date (t,-, 77,-), unde: ti => anul efectuării măsurătorilor rji => populația masurata in anul t Pentru aproximarea creșterii demografice este propus un model exponențial de forma: H(*i,*2) = xiexpX2t/, cu x = (xi,X2) reprezentând setul de parametrii ai modelului Notam r] = m]T si /W(x) = [/Wi(x) ■ ■ ■ /W,„(x)]r => estimarea parametrilor funcției creșterii demografice /W(x) se formulează ca o problema CMMP neliniara Frecvent, in aplicații de estimare si fitting, modelul jVt este liniar, i e de forma JVt(x) = Jx Jx = rj forma algebrica mjn xh- M2 x€Rn 2 forma optimizare In acest caz, funcția obiectiv f(x) = ^\\r] — _/x||2 este o funcție patratica convexa, avand Hessiana V2f(x) = JTJ 0 iar problema CMMP liniara devine: min f(x) ( = —1|?7 — _/x||2 xeR" v 7 \ 2 ' Cum f este convexa, atunci condițiile de optimalitate de ordinul I Vf(x) = 0 sunt si suficiente => orice soluție a sistemului JTJx — JTr] = 0 este punct de minim global a problemei CMMP De exemplu, daca rang(_/) = n, avem JTJ 0, deci: JTJx* -JTn = 0 x* = (jTJ^ 1 JTrj Definiție (Pseudo-Inversa Moore-Penrose) Fie matricea J G Rmxn de rang r, avand descompunerea valorilor singulare (DVS) data de J = WLV7 Atunci, pseudo-inversa Moore-Penrose J+ are expresia: J+ = VZ+UT, unde E = diag( n, II??—Jx||=minim forma algebrica v ■z forma optimizare Teorema In cazul general r = rang J media fj a punctelor rjj- Consideram setul de date {ti, , tm} si valorile corespunzătoare Dorim sa găsim vectorul parametrilor x = (xi,*2) a i- polinomul de ordinul intai p(t;x) = xi + x2t realizează predictia lui rj la momentul t Acesta poate fi găsit rezolvând problema: 1 X A 1 min, 5 12 = min, o xgR2 2 xgR2 2 i=l 1 tl unde J = _ 1 tn _ t2f? — tiȚt —tfj + rjt rj — J xi *2 2 Rezolvând problema se poate obține soluția , unde si at reprezintă variantele lui rj si resepctiv t, iar p coeficientul de corelație al acestora Consideram un sistem intrare-iesire: Pentru sistemul considerat dispunem de 40 de măsurători intrare-iesire y(t)}: Dorim sa aproximam sistemul printr-un model intrare-iesire de forma ARMA: ymodel(t) = Xlu(t)+x2u(t-l)+x3u(t-2)+X4u(t-3)+x5u(t-4) Gasirea setului de parametrii x = [xi ■ ■ *5]7 ai modelului poate fi realizata prin rezolvarea unei probleme CMMP liniare cu î? = [y(5) y(6) - y(40)]r: f u(5) "(4) "(3) "(2) "(1) \ "(6) "(5) "(4) "(3) "(2) J = u(7) "(6) "(5) "(4) "(3) \ "(40) u(39) u(38) "(37) "(36) / Răspunsul modelului estimat: Răspunsul estimat si răspunsul real odei (f) 0 5 10 15 20 25 30 35 40 In cazul JTJ inversabila, mulțimea soluțiilor optime X* conține un singur punct, i e X* = 1 Jrj^ Daca JTJ nu este inversabila, atunci o infinitate de soluții: X* = {xGR" : Vf(x) = 0} = {xGR" : JTJx - JTp = o} Ne interesează soluția din X* de norma minima => min illxll2 xeX* 2 Aceasta soluție este data de x* = J+rq si poate fi determinata rezolvând o problema regularizată CM MP liniara cu o constanta (3 > 0 suficient de mica: (=|lh M2 + (lM2 xER" \ Z Z Soluția problemei poate fi găsită prin rezolvarea sistemului obtinut din condițiile de optimalitate: Vf(x) = JTJx - JTr] + /3x= (jTJ + /?/„) x - JTr] = 0 ^x* = + J7/; Lema: Pentru o matrice J E Rmxn avem: lim (^J + ȘIn] 1 JT = J+ Demonstrație: Utilizând descompunerea DVS corespunzătoare matricei J = UHVT avem: (jrJ + /5/n) 1 JT = v(zTZ + /3ln^ 1ZTUT Evaluând partea dreapta a ecuației obținem: Ținând cont de definiția lui E+ si J+ se poate observa ca pentru —> 0 avem (JTJ + /?/„) 1 JT —> J+ Problema CMMP liniara poate fi interpretata in sensul determinării unui set de parametrii x gR" ce explica măsurătorile perturbate rj in cel mai bun mod Pentru 771, ,pm valori observate ale unei variabile aleatorii avand densitatea P(î?|x) ce depinde de setul de parametrii x, presupunem: î?/ = M,-(x) + unde: x => valoarea adevarata a parametrului x; fii => zgomot Gaussian de medie E(/3,-) = 0 si varianta E(AA) = ^2 fii si /3j sunt independente Introdusa de Fischer in 1912 gasirea estimatorului de verosimilitate maxima x* a lui x ce maximizează funcția de verosimilitate: m m P(jl\x) = 11 P(dlM = JJexp i=i i=i V 2^2 ) In general P(j]\x) si log P(j]\x) isi ating maximul in același punct x* => problema echivalenta: x* = arg max P(r/\x) = arg min — \og(P(ri\x}) xeR" xeR" = arg min |||S_1 (rj - M(x)) ||2, xCR" Z unde S = diag( ușor de rezolvat, R inf triungh Rezolvarea globala a problemelor CMMP neliniare in general NP-hard determinarea unui minim local se realizează iterativ Principiul de baza consta in aproximarea problemei la fiecare iterație cu liniarizarea sa in punctul curent Pentru probleme CMMP neliniare de forma: unde F(x) = rj — /W(x), metodele des utilizate sunt: metoda Gauss-Newton metoda Levenberg-Marquardt La fiecare iterație a metodei GN, fiind dat xk, se realizează următoarea actualizare: xk+1 = arg min -|| F(xk) + J(xk)(x - xk) ||2, xGRn z s s liniarizarea lui /-(x) in unde _/(x) = dF^ ■ Notând F(xk) = Fk, J(xk) = Jk si presupunând ca JkJk este inversabila avem echivalent: xk+i = xk + arg mm ^\\Fk + Jkd\\2 = xk - (jkJk^ JkFk Observație: direcția dk = — (dk Jk) 1 Jk Fk = —J^Fk este direcție de descreștere (Vf(x/t), dk) = —F^Jk(J^Jk)^1J^Fk oc avem dk —> 0 deci metoda nu aplica nici o corecție punctului curent In schimb, pentru fik —> 0, dk —> —J^Fk, deci coincide cu direcția din metoda G-N Metoda GN vs LM (k=4) Metoda GN vs LM (k=5) Metoda GN vs LM (k=6) Metoda GN vs LM (k=7) Putem observa ca gradientul funcției obiectiv f(x) = |||F(x)||2 a problemei CM MP neliniare este dat de următoarea relație: Vf(x) = J(x)rF(x) Pentru Vf(x) = 0 => direcțiile ambelor metode sunt nule, ceea ce este o condiție necesara de convergenta la un punct staționar Daca notam cu F,- componenta i a funcției F avem: m i=l Observam ca in cele doua metode, cel de-al doilea termen al funcției f este neglijat Daca funcția F(x) este aproape liniara sau F,-(x) sunt mici in apropiere de soluție atunci acest termen este redus De exemplu, daca se cauta o soluție a sistemului liniar F(x) = 0, cu m = n, acest termen este nul In plus, daca Jk este inversabila, atunci direcția Newton coincide cu cea din metoda Gauss-Newton Teorema: Fie x* un punct ce satisface condițiile suficiente de ordinul doi Iterația metodelor Gauss-Newton si Levenberg-Marquardt in apropierea lui x* este de forma = x/c — HfrVfțxk), unde Hk Q este data fie de Hk = (JTJk) \ fie de Hk = A + fîiJn) Presupunem satisfăcuta condiția Lipschitz: ||H4V2f(x/()-V2f(y)) || 7 > 0 Fiecare student alege o aplicație a i sa se modeleze ca o probleme de optimizare fara constrângeri Cerințe: descrieți in detaliu aplicația aleasa dati formularea matematica ca o problema de optimizare UNLP utilizați cel puțin 2 algorithmi de optimizare pentru UNLP dati la curs pentru rezolvarea ei dati codul Matlab pentru cei 2 algoritmi comparați rezultatele cu funcția Matlab fminunc comentați rezultatele obținute (grafice, tabele) Ion Necoara 2014 Problema de optimizare neliniara supusa la constrângeri: (NLP) : min f(x) xeR" s l : gi(x) „ ||x — x0^, i e X = : x>0 ([x°]>n)/ = 0 daca xf 0 V/ min II*- [1 1 ^lli xeR3 S l X1 + x2 + x3 3,/?>0 R s \ R + \\sj — x\\ 0 funcție obiectiv liniara constrângeri patratice convexe caz particular de (CP) Dispunem de n transmitatori, cu nivelurile de putere Pi, P2?•••? Pm si de n receptori fiecare receptor / recepționează semnalul transmitatorului / puterea receptata de receptorul / de la transmitatorul j este data de , unde Gjj reprezintă amplificarea pe linia de comunicație (/,j) Puterea semnalului primit de receptorul / de la transmitatorul / este • Zgomotul din rețea este dat de Smjn min Pi H \- Pn PeR" c I P ■ Smin Tomografie computerizata = tehnica ce folosește raze X (sau alte tipuri de radiații) pentru a produce imagini 2D/3D ale interiorului obiectului scanat Procedura de funcționare consta in: 1 Se achiziționează o serie de , din diferite unghiuri, ale obiectului scanat; 2 Prin intermediul proiecțiilor obținute, se reconstruiește interiorul obiectului cu ajutorul unui algoritm iterativ; www mathworks com In majoritatea cazurilor, radiațiile folosite sunt daunatoare; de aceea se urmărește achiziționarea unui număr minim de proiecții Formularea problemei: • Fie x G Rn imaginea interiorului de reconstruit • Pentru reconstrucție, dispunem de diferite măsurători liniare (proiecții) ale imaginii x: b; = A;x / = 1, ■ ■ ■ , m • Notam vectorul proiecțiilor b G Rm si A = [4-f ■ ■ ■ e Rmxn matricea de achiziție • Imaginea interiorului reprezintă soluția sistemului liniar (subdeterminat deoarece sunt mai puține măsurători m decât dimensiunea imaginii n): • Reformulare in termeni de problema CMMP: unde de obicei se alege a = 2 sau a = 0 sau a = 1 Alegem a = 0 V 1 pentru a induce o reprezentare rara a imaginii (vectorul soluție) Se dorește o reprezentare rara a imaginii, deoarece aceasta permite: compresie ușoara; algoritmi rapizi pt procesare; memorie de stocare mica; eliminarea ușoara a zgomotului; se da o matrice X cu elemente lipsa (e g o imagine) se presupune ca matricea X este de rang scăzut scopul este sa se gaseasca elementele lipsa din X pentru a impune rang mic asupra unei matrici se folosește nuclear norm || ■ ||*: daca A are descompunerea DVS A = ajUjvȚ, atunci ||4||* = ^-=1 07 problema se pune astfel: min rang(X) XgRmX" relaxare convexa min XgRmX" 11*11 s l : Xij = Ajj Vi,j G S s l : Xj = Ajj Vi,j G S unde se dau valorile Aj cu (/, j) G S o submultime a elementelor matricei cautate relaxarea convexa se poate scrie ca un SDP: min tr(Wi)+ tr(W2) x,wltw2 s L: Xij = Aj Vi,j G S, Wr XT X w2 >= o imaginea data (40% elemente lipsa) si imaginea recuperata Definim o măsură a raritatii ||x||o =» {i : x,- 0} si ||x||i = ^2,- |x,-| si apoi problemele de optimizare asociata sistemului +x = b: min llxllo s l Ax = b min llxlli s l Ax = b xeR" xeR" relaxare convexa la^a-l • Mutual coherence: u(A) = max ,, îmi u l Rp: min ||X||0 s l A(X) = b => min ||X||* s L X(X) = b X^KmXn X£RmXn relaxare convexa • Restricted Isometry Property de ordin r: cea mai mica constanta 5r{A) a i pt orice matrice X de rang cel mult r avem: 1 - 6r(A) 1 Atunci unica soluție a problemei originale neconvexe este matricea de rang cel mult r satisfacand A(X} = b • Presupunem ca r > 1 satisface d$r Rm, h(x) : Rn —> Rp Definim funcția Lagrange £ : Rn x Rm xRp ->R: Numita si Lagrangianul, funcția Lagrange a fost introdusa pentru prima data de Lagrange pentru probleme cu constrângeri de egalitate Variabilele: A G Rm, G Rp se numesc multiplicatori Lagrange (variabile duale) Regula generala: impunem ca multiplicatorii lagrange pentru constrângeri de inegalitate sa fie nenegativi, i e A > 0 Pentru orice variabila primala x fezabila pentru problema de optimizare (NLP) (i e g(x) 0, g(x) R este , pentru multiplicatorii A si p, fixați, i e: Daca q(A, /j,) este nemărginită inferior, i e q(A, /fi) = —oo, atunci spunem ca perechea (A, /z) este Lema: Pentru orice pereche duala (A,/t) fezabila, i e A G RȚ si fi G Rp, avem Demonstrație: pentru orice punct fezabil x avem q(Ă, fi) R este Demonstrație: observam ca Lagrangianul £(x, •, ) este funcție afina in (A, p,) pentru x fixat Fie a G , atunci pentru (Ai,^i) și (A2, ^2) avem: q(aAi + (1 — a)Ă2, aijUi + (1 — 0)^2) = inf £(x, aXi + (1 — a)Ă2, aijUi + (1 — 0)^2) xeR" = inf a£(x, Ai, ui) + (1 — a)£(x, A2, U2) xeR" >a inf £(x, Ai, ui) + (1 — a) inf £(x, A2,M2) xeR" xeR" = aq(Ai,^i) + (1 - a)q(A2, At2) Problema duala este definita ca fiind unde notam cu q* valoarea optima duala Observație: din moment ce q este intotdeauna concava (—q convexa), problema duala este intotdeauna problema convexa: q* = max q(X,u) = — min —q(X,u) A>0, ;teRP A>0, Observație: mulțimea fezabila duala Q = RȚ x Rp este o mulțime simpla (proiecția pe aceasta mulțime se poate realiza ușor) Fie următoarea problema patratica convexa: min l/2xrx xSR" s l: Ax = b, unde A G Rpxn Observam: f(x) = l/2xrx patratica convexa (Q = /„ G 0) si g(x) = 0 (fara inegalități) h(x) = Ax - b, h : Rn Rp, deci Lagrangianul £(x, ^4) = f(x) + //Th(x') este: =l/2xrx + pTțAx—b') => = min l/2xrx + p,T(Ax—b) \7xC(x, / 4) = 0 => x(^) + AT/i = 0 => x(^) = —AT11 qfai) = AAT / i — p,TAATp, — bT 11 = -1/2htAAth - bT/i Astfel, problema duala este de fapt (QP) fara constrângeri, dar Hesiana este pozitiv semidefinita! Fie următoarea problema CMMP cu restricții patratice: min 1/2|| Ax — b||2 xeR" s l: l/2||x||2 0) si /?(x) = 0 (fara egalitati) g(x) = 1/2||x||2 — c patratica convexa, deci Lagrangianul £(x, A) = f(x) + ^s(x) este: = 1/21| Ax — 6112 + l/2A(xrx — 2 c) => = min £(x, A) : Vx£(x, A) = 0 => ATAx(X) - ATb + A/„x(A) = 0 => x(A) = (AtA + Xln)-1ATb q(A) = £(x(A), A) = -l/2bTA(ATA + Xln)-1ATb - Xc + 1/21|£>||2 Astfel, problema duala devine (este concava!) Fie următoarea problema de optimizare NLP: 4 4 min X-, + x9 — 4xiX2 xeR2 s l: x > 0, xi + X2 = 5 Observam: f(x) = X| + x9 — 4xiX2, f : R2 —> R, n = 2 g(x) = —x, g : R2 —> R2, deci m = 2 si /?(x) = xi + X2 — 5, h : R2 —> R, deci p = 1 si Lagrangianul £(x, A, p) = f(x) + Arg(x) + pT h(x) este: £(x, A, p) = xf + X2 — 4xiX2 — AiXi — A2X2 + p(x-[ + X2 — 5) Astfel, problema duala este de fapt q* = max AeR^, ,ueR min xf + x9 — 4xiX2 — A1X1 — A2X2 + p(xi + X2 — 5) ,xeR2 > Din moment ce funcția duala este mărginită superior de funcția obiectiv, rezulta următoarea teorema: Teorema (dualitate slaba): următoarea inegalitate are loc pentru orice problema de optimizare (NLP): Observație: daca avem x* fezabil pentru primala si (A*,jU*) fezabil pentru duala astfel incat q(A*,jU*) = f* atunci x* este punct de de pentru primala iar (A*,jU*) este punct pentru duala Observație: daca f* = —oo atunci q(A,^) = oo,V (A, p,) G Q De asemenea, daca q* atunci problema primala este Diferența f* — q* se numește , iar din dualitate slaba este nenegativa In anumite cazuri (e g optimizare convexa in care mulțimea fezabila indeplineste anumite criterii), putem obține o versiune mai puternica a dualitatii Condiția Slater: presupunem ca problema (NLP) este convexa (i e f si gi, ■ ■ ■ ,gm sunt convexe, iar hi, ,hp sunt afine) si exista x gR" fezabil astfel incat Teorema (dualitate puternica): Daca problema (NLP) primala , atunci valorile optime pentru problemele primale si duale sunt egale, anume: Mai mult (complementaritatea), unde x* este punct de minim global pentru problema primala si (A*,jU*) este punct de maxim global pentru problema duala problema originala (primala) este echivalenta cu min ( max A, u) xeR" \A>0,/x problema duala este interschimbarea intre min si max max ( min C(x, X, u) ă>o,m \xeR" dualitatea slaba are loc intotdeauna deoarece maxmin£( ) 0 si 0 altfel (funcția indicator al lui R+) Definim lo(y) = oo daca y 0 si 0 altfel (funcția indicator al lui {0}) Rescriem problema originala (NLP) ca: min f(x) + 1-(^W) + $2 1o(/?j(x)) i j pentru a obține dualitate inlocuim l_(^(x)) cu A/^(x) (o măsură de disconfort când A,- > 0 si g/(x) > 0) De asemenea iijhj(x) este o margina inferioara pe lo(/?7(x)) mm+ 52 + 52 maw min(x — c — l)2 Funcția originala (albastru), constrângeri (verde), £(x, A) pentru diferiți A (punctat) s l 2 C Funcția duala oarea optima (punctat) q(A) si val-primala f* Pentru a vizualiza grafic diferența dintre dualitate slaba/puternica, consideram un caz particular de problema (NLP): minxe »{f(x) : g(x) - 0): f* = min -xTQx + qTx xgR" 2 s L: Cx — d 0, ;teRP 1 2 A A d+CQ-'q T A' b+AQ^q A |grQ_1g funcția duala este patratica in (A, /z) si concava problema duala este un (Hessiana nu mai este pozitiv definita) formularea duala a unei probleme QP are constrângeri mult mai simple: A > 0 si / i€Rp ultimul termen in funcția duala este o constanta, insa trebuie pastrat pentru a menține dualitatea puternica 0,xi — X2 > 1,xi + X2 = 3 Daca exprimam detaliat funcție patratica f(x) = |xrQx + qTx, f : R2 —> R, avem: rz A Q11 2 , Q22 2 , ^?12 + Q21 f(x) = —xx + ~^x2 d 2 -X1X2 + 91X1 + q2X2 Astfel, pentru problema noastra (dorim Q simetrica, Q12 = Q21), rezulta: Q = Q11 Q21 Q12 Q22 2 ii r 1' 1 2 ’ q ~ -1 Observam: A(Q) = {1,3}, Q >- 0, det(Q) = 3 si Observam: Q-1 1 3 2 -1 -1 2 g : R2 —> R2, g(x) = g(x) = Cx - d => C e R2x2, d & R2 ~*1 x2 - X1 + 1 -1 -1 0 1 d = 0 -1 h : R2 —> R, /?(x) = xi + x2 — 3 h(x) = Ax-b^Ae Rlx2, b e R, A = , b = 3 : se urmărește separarea unui set de obiecte in doua clase Modelare matematica: sa se determine (a) optim ce separa un set de vectori (obiectele) x-t in doua clase: clasa vectorilor x-, pt care aTXj — b —1 E-mail filtering - determinarea unui estimator ce separa un set dat de email-uri in doua clase: si dispunem de un dicționar de cuvinte D = {ci, C2, , cp}, unde ci reprezintă cuvântul i antrenam estimatorul printr-un set de e-mail-uri cunoscute pentru fiecare e-mail cunoscut i asociem vectorul xi = [«i, n2, • • •, nP]r, unde rij reprezintă numărul de apariții ale cuvântului cj in e-mail-ul i pentru fiecare e-mail cunoscut i asociem eticheta yi = {—1,1}, unde y,- = —1 daca e-mail-ul cunoscut i este spam, altfel y, = 1 min weR",beR,£>o s l yi(ivrxi - b) 0 si SLi ai = 1- Este clar ca cTv-, > f*, orice v,- fiind fezabil Notăm cu I mulțimea de indecși: T = {i : a,- > 0} Daca exista /'o G T astfel incat cTvi0 > f*, atunci: f* = cTx* = aiQcTViQ + a]CTVj > a/ f* = f* /ez\{/b} izz si deci obținem o contradicție Aceasta implica ca orice vârf pentru care at > 0 este un punct de minim Exprimam acum duala unui (LP) Lagrangianul asociata unui (LP) general este dat de: £(x, A, ju) = cTx + XT{Cx — d) + [iT {Ax — b) = —XTd — iJ,Tb+^c+CTX + ATij^ x Funcția duala corespunzătoare este: q(A,/z) —XTd — [jTb+ inf (c + CTX + ATfi\ x xeR" \ / —XTd — [iT b + 0 — CX3 dacă c + CTX + AT/j, = 0 altfel : problema din membrul drept de mai sus este un (LP) nemărginit, a cărui minim este — oc daca c + CTX + AT/j, 0 Astfel, putem exprima problema duala drept: 0, '-dl T FA' —b fi c + CTX + AT fi = 0 un LP poate fi scris in forma standard: min tcTx : Ax = b, x > 0}, xeR" prin folosirea de variabile suplimentare (numite si variabile artificiale) Duala sa rezulta: maxjtieR„{br/z : AT/i 0, X2 > 0 Astfel, problema vinarului poate fi formulata drept o problema (LP) cu constrângeri de inegalitate: min țcTx : Cx R, g : Rn —> Rm si h : Rn —> Rp sunt funcții diferentiabile de doua ori a problemei (NLP): X = {x G R" : g(x) 0,xi + X2 = 5 • Observam: funcția obiectiv: f(x) = x? + x^ _ 4xiX2, f : R2 —> R, n = 2 constrângeri de inegalitate g(x) = — x, g : R2 R2, m = 2 constrângeri de egalitate h(x) = xi + X2 — 5, h : R2 —> R, p = l • Mulțimea fezabila in acest caz X = {x G R2 : x > 0, xi + X2 = 5} este convexa! Teorema: Fie X o mulțime convexa si f G C1 Pentru problema de optimizare constrânsa min f(x) xeX avem următoarele condiții de optimalitate: daca x* este minim local atunci: Vf(x*)r(x-x*) > 0 Vx G X daca in plus f este funcție convexa, atunci x* este punct de minim daca si numai daca: Vf(x*)r(x-x*) > 0 Vx G X Punctele x* ce satisfac inegalitatea anterioara => pentru (NLP) avand constrângeri convexe Definiție: O constrângere de inegalitate g/(x) 0, xi + X2 = 5} mulțimea activa in punctul x = T este -4( r) = {!} Consideram intai probleme (NLP) ce au doar constrângeri de egalitate: (NLPe) : min f(x) xeR" s l : /)(x) = 0 Pentru problema (NLPe) toate constrângerile sunt considerate constrângeri active Vom studia condițiile de optimalitate pentru clasa de probleme (NLPe), pentru ca apoi sa extindem rezultatele obținute la clasa generala (NLP) O curba pe o suprafața S => o mulțime de puncte x(t) G S continuu parametrizate in t, pentru a mulțimea tuturor vectorilor x(t*) definiți de x(t) G S avand proprietatea ca x(t*) = x* Pentru o funcție h : Rn —> Rp, cu b(x) = [bi(x) bp(x)]r notam Jacobianul prin Vb(x) Reamintim ca Vb(x) este o matrice p x n ,avand elementul de pe poziția (/, j) egal cu Definiție: Un punct x* ce satisface constrângerea h(x*) = 0 se numește daca gradientii V/?i(x*), , V/?p(x*) sunt liniar independenți (Jocabianul are rangul p) Teorema: Intr-un punct regulat x* al suprafeței S definita de constrângerile de egalitate h(x) = 0, planul tangent este egal cu: M = {d G Rn : V/?(x*)d = 0} = {(/GR" : V/7,(x*)rd = 0 V/= 1, ,p} Lema: Fie x* un punct regulat al constrângerilor h(x) = 0 si punct de extrem local (minim sau maxim local) al problemei de optimizare (NLPe) Atunci, orice d G Rn ce satisface: V/?(x*)d = 0 trebuie sa satisfaca si: Vf(x*)rd = 0 • Fie constrângerea data de funcția h : R3 —> R, /)(x) = xf + xf + 3xi + 3x2 + X3 — 1 Consideram punctul x* = r pentru care h(x*) = 0 • Jacobianul lui /î(x) va fi: V/?(x) = iar V/?(x*) = are rang 1 => x* punct regulat • Din teorema precedenta pentru caracterizarea planului tangent, rezulta ca orice direcție tangenta d = [di c/2 ds\T va satisface V/?(x*)d = 0, deci 3 di + 3d? + d?, = 0 (un plan) Teorema: Fie x* un punct de extrem al funcției obiectiv f supusa la constrângerile /)(x) = 0, i e al problemei de optimizare (NLPe), si presupunem ca x* este un punct regulat pentru aceste constrângeri Atunci, exista un multiplicator Lagrange /j* G Rp a i : Punctele x* pentru care exista /_/ * a i condițiile (KKT-NLPe) sunt satisfăcute => pentru problema (NLPe) (puncte de minim, puncte de maxim sau puncte sa) Condițiile (KKT-NLPe) nu reprezintă altceva decât condițiile de optimalitate de ordinul I pentru minimizarea Lagrangianului neconstrans £(x, A) = f(x) + l^T h(x)- Vx£(x / i) = 0 si Vm£(x, ju) = 0 V£(x, A) = 0 Fie problema (NLPe): min x-y + X2 xeR2 s I : /z(x) = xf + x/ - 2 = 0 Putem observa ca orice punct fezabil este regulat Drept urmare, orice punct de minim local satisface sistemul V£(x,/z) = 0: 2/zxi = — 1 2/ZX2 = — 1 X1 + x2 = 2- Soluțiile (xj*,X2,jU*) ale sistemului sunt punctele staționare (—1, —1,1/2) si (1,1, —1/2) Sunt puncte de extrem cele doua soluții sau puncte sa Consideram următoarea problema de optimizare: min 4xi — 2x2 + xeR3 s l : /?i(x) = Xi — X2 + X3 = 0, /?2(x) = —xi + 2x2 — x| = 0 Jacobianul constrângerilor va fi de forma: V/?(x) = 1-11 -1 2 -2x3 ce are rangul minim 2 => orice punct fezabil este regulat Pentru a găsi punctele de extrem rezolvam sistemul dat de condițiile (KKT-NLPe): 4 + /zi — /z2 = 0, - 2 - /zi + 2/z2 = 0, 2x3 + Mi ~ 2x3/z2 = 0 xi - x2 + x3 = 0 —xi + 2x2 — *3 = 0 Obținem ca soluția unica a sistemului este data de x* = [—1 0 1] si Mi = —6, M2 = —2 Care este natura punctului x*? Consideram problema: min — x-y xeR2 s l : /7i(x) = (1 - xi)3 + x2 = 0, /72(x) = (1 - xi)3 - x2 = 0 Putem observa ca problema are un singur punct fezabil x* = r => minim global De asemenea, avem Vf(x*) = [—1 0]T, \7hy(x*) = r si V/?2(x*) = r => x* este un punct regulat! In acest caz, condițiile de optimalitate de ordinul I nu pot fi satisfăcute => ^/zi si /z2 a i : 0 1 /zi + 1^2 0 -1 1 o Teorema: Presupunem ca x* este un punct de minim local al problemei (NLPe) si un punct regulat pentru constrângerile aferente Atunci exista /z* G Rp astfel incat Vf(x*) + V/?(x*)V = 0 si h(x*) = 0 In plus, daca notam planul tangent in x* prin M = {d E : V/7(x*)d = 0}, atunci matricea Hessiana a Lagrangianului in raport cu x p V2£(x*,/z*) = V2f(x*) + J2/z*V2/?/(x*) i=i este pozitiv semidefinita pe planul tangent M, adica drV2£(x*,/z*)d > 0 VdG/W Consideram din nou problema: min x-y + X2, xgR2: /7(x)=x2+x22-2=0 pentru care matricea Hessiana a Lagrangianului in variabila x este: r ? o V2£(x, [1) = V2f(x) + /zV2/?(x) = fi 0 2 • Ținând cont ca planul tangent la un punct x^O este de forma M = {d : d = 0}, pentru care putem alege o baza D(x) = [~x2 xi]T avem: D(x)rV2£(x,//)D(x) = 2//(xf + x2) Pentru prima soluție a sistemului V£(x,/j) = 0 avem d1rV2C(—1, —1, l/2)di = 2 > 0, unde dy = D(—1, —1) => posibil punct de minim local! Pentru cea de-a doua soluție d2rV2£(l, 1, — l/2)c acesata soluție nu este minim local! Teorema: Presupunem un punct x* G Rn si un [i* G Rp a i : Vf(x*) + V/?(x*) V = 0 si h(x*) = 0 Presupunem de asemenea ca matricea data de Hessiana Lagrangianului V2£(x*, //*) = V2f(x*) + Z1/V2/?,(x*) este pe planul tangent M = {d : V/?(x*)d = 0}, adica drv2z:(x*,/z*)c/> 0 VdG/W\0 Atunci x* este punct de minim local strict al problemei (NLPe) Fie problema: min —X1X9 — *i *3 — *2*3, xeR3: xi+x2+x3=3 pentru care condițiile de optimalitate de ordinul I conduc la următorul sistem liniar: -x2 -x3+ fi = 0 -*i -x3 + fi = 0 -x3 - x2 + fi = 0 x3 +x2+ x3 = 3, avand soluția unica x^ = x^ = x^ = 1 s\ fi* = 2 Hessiana Lagrangianului va avea următoarea forma: V|£(x //) = V2f(x) = 0 -1 -1 -1 0 -1 -1 -1 0 Planul tangent intr-un punct x fezabil va fi de forma M = {d : d = 0} => di = —c/2 — ds => alegand [di d?] = si [di d?] = [—1 — 1] => o baza a planului tangent la suprafața definita de constrângerea d(x) = 0 va fi: D(x) = 0 2 1 -1 -1 -1 Drept urmare, obținem ca x* este minim strict local intrucat: D(x*)rV|£(x* /r)D(x*) = 2 0 0 2 >- 0 Este Hessiana funcției obiectiv evaluata in x* pozitiv definita? In continuare, vom extinde condițiile de optimalitate derivate anterior la cazul general: (NLP) : min f(x) xeR" s l : g(x) O si g(x*) 0 implica g/(x*) = 0 Mai departe, ținând cont ca x* este un punct de minim local pentru problema (NLP) Notând cu X mulțimea constrângerilor problemei (NLP), atunci x* este un punct de minim local si pentru problema avand mulțimea de constrângeri o submultime a lui X si definita prin setarea constrângerilor active la zero Astfel, pentru problema (NLPe) ce ar rezulta, definita pentru o vecinătate a lui x*, exista multiplicatori Lagrange si deci prima relație este satisfăcuta pentru X* = 0 daca g/(x*) 0 si drept urmare si cea de-a doua relație este satisfăcuta Mai trebuie aratat ca A* > O pentru constrângerile active g/(x*) = 0 Fie o componenta X*k 0 suficient de mic x(t) este fezabil iar din prima relație a condițiilor (KKT) obținem: ——— t=o = Vr(x )d „ £(x, A*,jU*), reprezentând condițiile de optimalitate de ordinul I, adica Vx£(x*,A*,^*) = 0 Cea de-a doua relație reprezintă condițiile de : g(x*)rA* = 0 g/(x*)A* = 0 V/' Ultimele doua relații exprima fezabilitatea primala (g(x*) 0) Consideram problema: min 0 5 (x? + x? + x?) xeR3 v ' s l gi(x) = Xi + x2 + x3 + 3 X1* + x2 + x3 Aî = A2 = 0 => x/ = x2 = x3 = 0 => contradicție cu ipoteza; 2 => x* + x* + x* 0 => Ai = —A2 = 0 => x2 = x3 = 0 => contradicție cu ipoteza; 3 => x/ + x2 + x3 = —3, x/ 0 si Ă2 = 0 => putem alege x/ = x2 = x3 = —1 si A3 = 1 ce satisfac condițiile (KKT) => ; 4 => x3 + x* + x* = -3, Xi = 0 si Aî, A2 > 0 => x2 = x3 = —3/2 si Ai = —A2 = 3/2 contrazice condiția A2 > 0 Consideram problema: min 2xf + 2xiX2 + xf — 10xi — 10x2 xeR2 s l : xf + xf 0 Pentru acest exemplu putem considera doua constrângeri active, una sau niciuna Presupunem prima constrângere activa iar cea de-a doua inactiva, rezultând un sistem de trei ecuații corespunzător, avand soluția xf = l,xf = 2, Af = 1 si Af = 0 Soluția verifica 3xi + x2 1 => satisface condițiile (KKT) Planul tangent in punctul regulat x* pentru problema generala (NLP) este planul tangent corespunzător pentru constrângerile active: M = {d : Vgj(x*')Td = QVj G A(x*), Vhitx^d = 0 V/ = 1, , p} Teorema: Fie f, g si h funcții continuu diferentiabile de doua ori si x* un punct regulat pentru constrângerile din problema (NLP) generala Daca x* este un punct de minim local pentru problema (NLP), atunci exista A* G Rm si p* G Rp a i condițiile (KKT) sunt satisfăcute si in plus Hessiana Lagrangianului in raport cu x: p m V2£(x*, A*,p*) = V2f(x*) + J>*V2/7/(x*) + A*V2g,(x*) i=i i=i este pozitiv semidefinita pe subspatiul tangent al constrângerilor active in x*, adica: drV2£(x*, A*,p*)d > 0 Vd £ M Teorema: Fie f,g si h funcții continuu diferentiabile de doua ori Fie de asemenea: un punct regulat x* G Rn variabilele duale A* G Rm si /z* G Rp pentru care condițiile (KKT) sunt satisfăcute nu avem constrângeri de inegalitate degenerate, adica AJ > 0 pentru orice j G Al(x*) Hessiana Lagrangianului V2/}(x*, A*,/z*) este pozitiv definita pe planul tangent M Atunci x* este un punct de minim local strict pentru problema (NLP) generala Fie problema: min X2 xgR2: x2+x2-l 0 => constrângerea este activa Adaugand aceasta condiție la cele doua anterioare, obținem soluția x* = r si A* = 1/2 Planul tangent in x* este dat de {d : d = 0} = {d : d? = 0} De asemenea, Hessiana Lagrangianului V2£(x*,A*) = 2A*/2 este pozitiv definita x* punct de minim strict Teorema: Fie o problema convexa (CP) de forma: (CP) : f* = min f(x) xeR" s l : g(x) x* = AT(AAT)~1 b = A+b, unde A+ este pseudoinversa lui A Extindem problema anterioara la cazul funcțiilor patratice generale: 1 t T min -x Qx + q x, x: Ax=b 2 unde Q >- 0 Condițiile (KKT-CP): Qx* + q + AT[j* = 0, Ax* = b Obținem astfel ca soluțiile primale si duale au următoarele forme: /i* = -(AQ~1AT)~1[AQ~1q + b] x* = -Q^A7/!* - Q~xq Consideram următoarea problema convexa: min (xi — 5)2 + (x2 — 5)2 xgR2 s l : gi(x) = xf + xf - 5 = A4 = 0 => condițiile (KKT-CP) se reduc la următorul sistem: 2(x*-5) 1 r 2x* n r aî 2(x* - 5) 2x* 1 H Ă2 (xf)2 + (x*)2 - 5 = 0, X* + 2x* -4 = 0 Obținem soluția x* = r => x* punct de minim global Ion Necoara 2014 Problema de optimizare neliniara supusa la constrângeri (A// P) : f* = min f(x) xeX f diferentiabila (nu neaparat convexa), X mulțime nevida, inchisa si convexa Exemple: min xeR" s l I d Reamintim condițiile de optimalitate de ordinul I: min f(x) xSR" Sub presupunerea ca f G C1, orice punct de minim local satisface sistemul: min f(x) xeX Sub presupunerea ca f G C1 si X mulțime convexa, orice punct de minim local x* satisface inegalitatea: Demonstrație: Presupunem ca exista y G X astfel incat Vf(x*)r(y-x*) 0 exista 0 G astfel incat: f(x* + t(y - /)) = f(x*) + tVf(x* + 0t(y - x*))r(y - x*) Din continuitatea lui Vf, alegem t 0 si x E X Observam ca d = ^2 0 & xk G X deci xk+i = xk + ak(xk ~ xk) Clar, daca Xk nu este staționar, atunci exista Xk G X a i — Xk) 0 Teorema Fie un sir Xk generat de metoda direcțiilor fezabile = Xk + ctkdk- Presupunem ca direcțiile dk sunt conectate prin gradient la Xk pasul (%k este ales prin metoda ideala sau backtracking Atunci orice punct limita al șirului Xk este Demonstrație: Consideram alegerea lui ctk prin backtracking (cazul ideal este similar) Presupunem ca Xk converge la un punct nestationar x Șirul f(xk) este descrescător si convergent la f(x), i e f(xfc) - 0 Condiția de backtracking: f(xk) — f(x/(+i) > —cia/fVf(x/f)rd/f, implica akVf(xk)Tdk —> 0 Direcțiile dk sunt conectate prin gradient la Xk, deci dk este mărginit si lim sup Vf(x/f)rd/f o ceea ce implica: Din aceasta relație si procedura backtracking avem ca exista k a i f(*k) ~ f(*k + (afc/p)4) k Din teorema valorii medii avem ca exista ctk G a i -Vf(xk + ăkdk)Tdk k Considerând ca dk este mărginit (are un subsir convergent), evaluam limita in relația precedenta si obținem 0 R diferentiabila min xex Metoda Gradient Proiectat: Interpretare */ 0, Vx G X s__________ak______________, Vq(xfe+i) Considerând x = xk si prelucrând condiția, avem: Vf(x/()r(x/(+i -xk) 0 sunt daca X compacta si deci putem folosim rezultatul precedent de convergenta min f(x) = xeR2 V ? ! x + [-l — 1] x S l 2xi + X2 = 1 Observații: V f (x) = Qx + q = 2xi + x2 - 1 xi + 2x2 - 1 proiecția pe hiperplan: 2 1 Fie y G R2, [y]x = y - (2^1+^-1) Metoda Gradient Proiectat: x/j+i = [x/j — akțQxk + q)]x Pornind din xq = r si considerând ao = 1 avem *1 = [*o - <>o(Qxo + q)]x = 3/5 ' -1/5 miiiimize — — b | s t ||a; 1 IR diferentiabila si mulțimea X convexa si compacta (inchisa si mărginită) minf(x) xGX obiectiv liniara; Recomandata in cazul in care complexitatea per iterație este mult mai scăzută decât cea necesara pentru rezolvarea problemei originale Interpretare xk = argmin f(xk) + Vf (xk)T(x - xk) x£ X adica folosește doar o aproximare liniara pentru funcția f in comparație cu metoda gradient proiectat care folosește o aproximare patratica pentru f cu Hessiana ^ln Gradientul condițional este caz particular de metoda fezabila: ^f(xk)T(xk -xk) R diferentiabila si X mulțime convexa min f(x) xeX Introducem următoarea norma vectoriala: H >- 0 => ||x||h = y/xTHx Metoda Newton Proiectat: norma indusa de o matrice H unde presupunem ca Hessiana V2f(x/ - 0 (pozitiv definita)! Interpretare xk+1 = argmin 0, Vx G X s__________&k_____________________, Vq(xfe+i) Considerând x = Xk si prelucrând condiția avem: Vf(xfe)rQxfe+1 -xfe) 0 sunt daca X compacta si folosim rezultatul precedent de convergenta min xgR2 f(*) = | IMI2 - log(xi ~ ^2 + 1) cTx-\-d S l *1 + X2 = O aTx=b Metoda Newton Proiectat: xk+1= arg min V f(xk)T(y - xk) + ^-(y - xk)T^2f(xkXy - xk) aTy=b K Observații: W) = x - = * - ^b+T V2f(x) = /2 + {cT,+dyCCT = I 1 [1 1 (X1-X2+1)2 -1 1 Soluție explicita a iterației Newton (cu pas ak = 1) in doua moduri: 1 Utilizând dualitatea Xk+i-Xk (v aT^2f^-la(y2f(Xk^ la)’ unde v = (V2f(x/()) 1Vf(x/()- 2 Eliminând constrângerea liniara, i e yi = y- \ b — ajYj 1 \ J xk+1 = arg min Vf(x/ 0 astfel incat daca ||xo — x*|| R si h : Rn —> Rp sunt funcții de doua ori diferentiabile Condițiile de optimalitate de ordinul I pentru aceasta problema (condițiile KKT): fie x* punct de minim atunci exista /z* G Rp astfel incat condițiile (KKT-NLPe) au loc: V£(x*,/z*) = vxc(x*,/z*) VMC(x*,/z*) o unde £(x,/z) = f(x) + pTh(x) este Lagrangianul Reamintim ca x* trebuie sa fie punct regulat, i e rangul lui V/?(%*) este p Fie următoarea problema de optimizare: min 4x,2 + 2x| + 4xiX2 — *1 + *2 xeR2 s l: xi — X2 = 1 f(x) = 4x2 + 2x| + 4xiX2 — xi + X2, f : R2 —> R /)(x) = xi — X2 — 1, h : R2 —> R £(x, = 4x2 + 2xj + 4xix2 - xi + x2 + /z(xi - x2 - 1) Din condițiile de optimalitate vxr(x*)Ă?) VM£(x*,/z*) = 0 rezulta sistemul: 8x1 + 4x2 — 1 + // = 0 4x2 + 4xi + 1 — // = 0 Xi — X2 — 1 = 0 Consideram sistemul dinamic liniar discret: zt+i = Atzt + Btut Vt > 0, unde zt G F?’z este starea, si ut G Rn“ intrarea sistemului Consideram costuri de etapa: £ț(zt) = l/2||zt — Zțef\\^t, £UM = 1/2|| ut-uf\\2Rt Formulam problema de control optimal N N-l min]Tvf(zt) + Zt,ut t=l t=0 s l : zt+i = Atzt + Btut Vt = 0, , N — 1 Definim variabila de optimizare: x = [uor zf ul uTn_i zl\T G R^+^) si vectorul referința aferent x”f = [Wf)T WT • • «1/ (4*f)T Atunci funcția obiectiv poate fi scrisa drept: f(x) = ^xTQx — qTx unde Q = diag(/?0, Qi, • • •, Rn-i, Qn), iar q = Qxref problema este de fapt un QP avand constrângeri de egalitate Constrângerile din problema de control optimal sunt constrângerile de egalitate ce provin din dinamici: zt+i = Atzt + Btut Prima constrângere poate fi scrisa drept: — BqUq + lnzZi = AqZq A doua constrângere poate fi scrisa drept: —A-[Zi + Bi Ui + lnzZ2 = 0 Considerând astfel dinamicile peste intreg orizontul de predictie, atunci constrângerile pot fi scrise drept: /)(x) = Ax — b = 0 unde -Bo 0 0 0 0 0’ AqZq 0 -Al -Bi Inz •• 0 0 0 si b = 0 _ 0 0 0 0 • -An-i -Bn-i lnz 0 Forma standard problema QP cu constrângeri de egalitate: (QPe) : min ^-xTQx — qTx v ’ xeR" 2 s l : Ax = b, Daca exprimam condițiile de optimalitate V£(x*,jU*) = 0, rezulta Qx — q + At,i = 0 Ax = b sau in notatie matriceala (numim K Q AT] Fx’ A Q q = Ky = q b Observație: matricea (KKT) este Daca matricea A G Rpxn are si pentru orice d G kernel(4) cu d 0 avem , atunci matricea (KKT) este (demonstrații) Daca A are rang p, atunci orice punct xGR" este Mai departe, daca pentru orice d G kernel(A) cu d 0 avem dTQd > 0, atunci condițiile suficiente de ordinul II sunt satisfăcute pentru (QPe) In concluzie, pentru o problema patratica cu constrângeri de egalitate în forma (QPe), existenta unui minim local este echivalenta cu inversabilitatea matricei (KKT) Rezolvarea sistemului KKT, Ky = n, se poate face prin următoarele abordări: presupunem Q inversabila si A are rangul p, : x = —Q^țĂTp, — q), introducem AQ_1 ATp = AQ_1 q — b p —> * Observam ca este inversabila, daca A are rangul p factorizarea LU K => utilizează o factorizare Cholesky : PTKP = LDLT, unde P este o matrice de permutare, D este o matrice bloc diagonala (dim bloc: 1 sau 2) si L este inferior triunghiulara metoda spațiului nul gaseste o baza (coloane ale lui) Z G Rnx(n-p) pentru kernel(4) si definește x = Zv + y, unde y este o soluție particulara a constrângerilor Ay = b putem elimina constrângerile Ax = b reformulam problema originala ca o problema QP (x —> Zv + y): min |(Zv + y)T Q(Zv + y) - qT(Zv + y) |/SR"-P 2 Daca sistemul KKT nu se poate rezolva explicit precum in cazul problemelor QP, atunci trecem la metode iterative pentru NLP-uri cu constrângeri de egalitate Metodele Lagrange se bazeaza pe condițiile KKT: Vx£(x, f) = 0, Kx) = 0 • V£(x, f) = 0- Definim: F(y) = V£(x,//) = Soluțiile problemei NLPe se sistemului neliniar F(y) = 0 Forma generica a metodelor Lagrange: Iterațiile metodelor Lagrange depind Pasul ctk poate fi stabilit conform unei Observație: 7-"(x, ^)>0, J7(x //) = 0 Vx£(x, a) = 0 si /?(x) = 0 Astfel, daca rezolvam problema neconstransa: min 7(x,u) xeR",/-teRp (= |l|vxr(x,(J)||2 + jmil2) (i) => orice (x*,^*) al acestei probleme V£(x*,^*) = 0 Lema: presupunem ca (x*,^*) este punct de minim local al problemei (1) care satisface condițiile necesare de ordinul I Presupunem de asemenea ca rangul lui V/?(x*) este p si Hessiana Vx£(x*,jU*) este pozitiv definita Atunci (x*,/i*) este punct de minim global pentru (1), adica = 0 FM = V£(x,p) = vx£(x, At) /,(x) = o Metoda gradient pentru rezolvarea sistemului de ecuații F(y) = 0 O yk+1 = yk — aF(yk) Daca yk converge la y* atunci F(y*) = 0 Metoda Lagrange de ordin I are iterația: Iterația se poate scrie : y/j+i = y/ - 0, atunci — 1^Lklk 0 Iterația (ML-I) va converge către un punct ce satisface VX£(x*, /j,*) = 0 dar nu putem garanta ca h(x*) = 0 Putem imbunatati convergenta (ML-I) prin alegerea altei funcții merit si pasul ak corespunzător acesteia: = |||Vx£(x,/z)||2 + |||/7(x)||2 -7£(x,/z), cu 7 suficient de mic Sub ipoteza Lk = V2x£(xk,iik) 7 0, (-Vx^(x/c, jUk), /j(xfc)) este o direcție de descreștere pentru funcția merit Jy [{Lk!k + Ăkhk-^lk)T (Aklk -^fhk)T][-lk hk]T = ~lk (Lk - ^ln)lk - ~fh2k Rn si H : Rnxp —> Rp funcții diferentiabile astfel încât x* = L(x*,/i*) si p* = H(x*,p*) Mai mult, presupunem ca toate valorile proprii ale matricei de dimensiune (n + p) x (n + p) VxL(x*,/d*) VxH(x*,p*} sunt în interiorul cercului unitate Atunci, (x*,p*) este un punct de atracție al iterației (ML) si când șirul generat {xk, p,k) converge la (x*,p*), rata de convergenta este liniara Avem astfel următoarea teorema pentru metoda Lagrange de ordinul I: Teorema: presupunem ca f si h sunt funcții de doua ori diferentiabile, x* este punct de minim local pentru care exista /z* satisfăcând condițiile (KKT-NLPe) Presupunem de asemenea ca x* este regulat si Hessiana V|£(x* //*) este pozitiv definita Atunci există ă > 0 astfel încât pentru orice a G (0 ă], (x*,/z*) este un punct de atracție al iterației (ML-I) si daca șirul generat (x/q/z/ sistemul are soluție Obținem următoarea iterație Newton clasica: unde direcțiile (dk,d^) sunt ~V2x£(xk,vk) (Nh(xk))T~ dk Vx£-(xk, f-lk) Vh(xk) 0 -h^xk) Rezultatele standard ale metodei Newton sunt aplicabile si aici: daca punctul inițial (%o,^o) este , atunci șirul generat de metoda Newton (ML-N) este convergent si converge către soluție cu rata Fie următoarea problema de optimizare: ■ 4 4 9 9 min x-i + x9 — x-i x9 T xi T x9 xeR2 s l: xi + x9 = 2 Vx£(x/(,/Ifc) — V2£(xk,Vk) = 4x^ — 2x^X2 + 1 + /j 4X2 ~ ^xlx2 + 1 + // ’ 12x| — 2x| —4x^X2 —4X1X2 12x2 ~ ^X1 : in general avem h : Rn —> Rm, iar V/?(x) reprezintă funcției In cazul nostru, avem /î(x) = xi + x9 — 2 iar Jacobianul rezulta Vh(x) = Rezulta astfel matricea KKT: 12x2 — 2x^ —4X1X2 1 —4xiX2 12x9 — 2xi 1 1 1 0 Avem astfel sistemul: 12x/ — 2x2 —4X1X2 1 —4X1X2 12x2 — 2x/ 1 dk ^k -Vx£(xk,Hk) -h(xk) 1 1 o Consideram un punct inițial (xo,/zo) cu xo = r si /zo = 0 înlocuind in sistemul anterior, rezulta: Rezulta: 12 0 1 -0 —2 1 1 1 0 do kJ -5 -1 1 do kJ —3/5’ 8/5 11/5 F2/5! (xi //i) = (x0 //0)+ (d0 d/) = 8/5 11/5 Reamintim notațiile: Lk = Vx£{xk, /J,k), lk = \?x£{xk, p,k), hk = h(xk) si Ak = Vh(xk) Astfel, produsul scalar dintre direcțiile Newton {dk, si gradientul funcției merit J7 este (ținând cont de sistemul KKT precedent satisfăcut de direcțiile {dk d^))'- [Lklk+A^hk Aklk][dk {dk")T]T — lk ^-kdk + hiAkdk + lk Ak d{ = -HM2 - IIMI2- Aceasta expresie este (adica condițiile (KKT-NLPe)) Astfel, direcțiile Newton {dk,d^) sunt Metoda Lagrange-Newton este o metoda de descreștere pentru o funcție merit —> are proprietăți de convergenta globala când iterația se ia cu un pas variabil Definim astfel metoda Lagrange-Newton generala prin următoarea iterație: Pasul ak este , i e ak = arg min F(xk + adk,^k + atf) Daca adaugam expresia (V/?(x/ R funcție convexa si A G Rpxn are rangul p Condițiile (KKT) pentru aceasta problema sunt , i e x* este punct de minim daca si numai daca exista p* G Rp astfel încât: (KKT - CPe) : Vf(x*) + AT p* = 0, Ax* = b In problema (CPe) putem elimina constrângerile de egalitate pe baza observației: {x : Ax = b} = {Zv + y : v G Rn~p{, unde coloanele lui Z reprezintă o baza pentru kernel(4) si y este o soluție particulara a sistemului Ax = b In acest caz putem rezolva problema de optimizare convexa fara constrângeri: min f(v) (=f(Zv + y)) i/SR"-p Daca v* este o soluție a problemei convexe fara constrângeri, atunci x* = Zv* +y este o soluție a problemei originale Fie Xk fezabil, următorul punct al iterației Newton se calculează folosind direcția Newton dk ce rezulta din aproximarea patratică de ordinul II a lui f în jurul lui Xk, i e c4=argmin f(d) (= f(xk) + Vf(xk)Td + ^dTV2f(xk)d') s l : A(xk + d) = b (2) Aceasta problema este un QP convex având constrângeri de egalitate Ad = 0, —> direcția Newton este caracterizata de sistemul liniar: V2f(xk) AT~ A 0 dk Pk+l -vfM -o Introducem decrementul Newton: 1 /o i'(xfc) = ^d^2f(xk)d^ : — f(dk) = 'y(x/ 0 Consideram L cantitatea (in bytes) a datelor (A,b,c) • Metoda simplex dezvoltata de 1947-1951; performante practice bune; complexitate in (n, / ) - Klee, Minty (1972) Problema: Exista algoritmi polinomiali (in n si / ) pentru probleme LP? • Metoda elipsoidului (punct interior) dezvoltata de 1979; performante practice slabe; complexitate per iterație: O(n2/_); complexitate : O(n4/_); • Metoda de punct interior dezvoltata de 1984; performante practice superioare metodei elipsoidului; complexitate per iterație: O(n/_); complexitate : O(n3-5/ ); Nansndm Karmarkar Problema neliniara convexa: min xSR" f(x) s t Ax = b, g(x) (y/î7 ln(z^/e)); pachete software foarte bune: CPLEX, MOSEK, CVX, MATLAB metoda simplex pentru LP are performanta practica buna dar in worst case complexitate exponențiala De ce? In medie, numărul de pivoți este: O(n2) netoda de punct interior pentru LP: in medie numărul de iterații este O(y/n In n) i e cel mult 30 — 50 de iterații pentru a găsi o soluție (NLP) : min xeR" s l : unde f : Rn —> R, g : Rn —> Rm si doua ori diferentiabile h : Rn —> Rp sunt funcții de min xeR" s l I d Condițiile de optimalitate de ordinul I (condițiile KKT): Vf(x) + V/?(x)r/z + Vg(x)rA = 0 g(x) 0 /7(x) = 0, unde A = diag(Xi, ■ ■ ■ , Am) Rezolvare dificila prin metode uzuale (e g metode mulțimilor active): condițiile de complementaritate Se aproximează cu un set de condiții netede Condiții de optimalitate pertubate (condiții KKT-IP): Vf(x) + V/?(x)rAt + Vg(x)A = 0 g(x) 0 /7(x) = 0 Pentru t = 0 obținem sistemul KKT original Algoritmii de punct interior rezolva succesiv sistemul (KKT-IP) pentru o serie descrescătoare de valori > 0 Soluția sistemului (KKT-IP) (x(r), A(t), ^(r)) converge la soluția sistemului (KKT) (x*,A*,jU*) când r —> 0 (cale centrala) Problema (NLP) se rescrie echivalent: min f(x) + xgR" s l /)(x) = 0, unde /_ : R —> R, /_(x) = daca x 0 1T 1T Funcția indicator /_(■) este nediferentiabila; Căutăm aproximare cu funcția bariera logaritmica; min f(x) — xgR" s t /)(x) = 0 = |x| - 1 /-(g-(x)) = f°> daca |x| daca |x| g(x) = x L(g(x)) 0 daca x daca x Condiții de optimalitate ale formulării cu funcția bariera: m -12 w+v/? wT/i =0 /=1 v— A; /j(x) = 0 1T 1T Avem g,(x) 0 obținem condițiile (KKT-IP) m vf(x) + E A'W/« + v/7(x)rĂ( = o i=l g(x) 0 /7(x) = 0 min xeR" f(x) s l g(x) Notam funcția bariera £>(x) = — log(~&(x)) s' punctul de i=i optim x(t) ce satisface si => Punctele de optim x(t) definesc si satisfac Vf(x(r)) + VB(x(t)) + AT/i = 0 => Orice punct x(t) produce un punct dual fezabil Alegem punctul dual fezabil definit de A/(r) = s' ju(r) = rp, Funcția duala in acest punct satisface m t7(Ă(r), ju(t)) = f(x(r)) + A,(T)g,(x(r)) + li(t)t(Ax(t) - b) i=i = f(x(r)) — mr O, a 0 Cat timp rriTk > e repeta: Calculează Xk+i = x(jk) pornind din punctul inițial Xk ("warrn start”); Descrește parametrul Tk+i = => Pentru a determina x/t+i rezolvam (e g prin Metoda Newton) reformularea cu bariera logaritmica cu parametrul Tk, pornind din punctul inițial Xk; => După k iterații avem f(xk) — f* (x) = — ȘA log(bj — A;x) prin x(t) = {x : B(x) = B(x(r))} mulțimile nivel C,T aproximează din ce in ce mai bine frontiera mulțimii fezabile {x : Ax 0 LP cu m = 100 inegalități si n = 50 variabile du^lity gap pornim cu x de pe calea centrala (to = 1) terminam când r = IO-8 pentru calcularea centrilor x(t) utilizam metoda Newton pentru egalitati folosind backtracking numărul total de iterații Newton pentru un a fixat este de ordinul zecilor numărul total de iterații Newton nu este foarte senzitiv pentru a 0) in cost utilizând funcția bariera logaritmica: min xSR" f(x) s l g(x) = s, h(x) = 0 (KKT-IPs): Vf(x) + V/?(x)rĂ( + Vg(x)A = O g(x) + s = O, h(x) = O s > O, A > 0 unde A = diagA => MPI rezolva sistemul perturbat prin metoda Newton; Fie x punctul curent, obținem Vf + ^hT/i + \7gTX /\s — re h g + s, unde S = diag(s) si £(x, s, A, /i) = f(x) + (g(x) + s)rA + h(x)T/i Iterație Metoda Newton: xk+l — xk + (tk^k ■> P-k+1 = P'k + ak^k’ $k+l — $k + pasul ak se alege in intervalul (0, «{fax]; => pasul axk se alege pe baza unei m m = f(x)-r^log(s/)+z/||/7(x)||i+z/^|max{0, g7(x)}| i=i i=i Se aleg p, q G (0,1], pasul axk = pm, unde m este cel mai mic intreg ce satisface relația (de backtracking) M„(xk + pmdk, sk + pmdk, Tk) - 0 si x, y vectori cu dim compatibile: IMI1/ = xTUx xT Ux + yT Vy = [xT yT] °1 U 0 • Reformulam funcția obiectiv utilizând observația precedenta: 1 Fii ,,ref 1 F P 1 U/-1 — Ui_1 K,_i 2 z-,- zf 0 /=1 L / J L 01 Q/J [ Zi-Z^ = |(x - xrefy Q(x - xref), unde am folosit notațiile: X=[uorz1r ■■■UT_1ZT], si Q = diag(Ro, Qi, ■ ■ ■ , Rn-i, Qn)- Formularea problemei de control optimal pe orizont finit rara (fara eliminarea stărilor): min S l 1T 1T unde forma standard de programare patratica matricile QP-ului sunt rare: Q si C sunt bloc diagonale, iar A bloc tridiagonala se aplica algoritmi eficienți specifici optimizării constrânse in funcție de dimensiunea problemei Dorim reformularea dinamicilor prin eliminarea stărilor, i e variabila devine: x = G zi = 4zz0 + Buuo Z2 = >4zzi -f- Buui = A^zq + AzBuuq + Buui ZN — Azz^—i + Buun—i — A^zq + A? 1 Bliuq + ■ ■ ■ + Buun—i => Notând z = [zf ■ ■ ■ z^]r G RWziz obținem: Obținem forma QP, reformuland constrângerile pe stare/intrare => Concatenam constrângerile pe stare: CzZi Combinând cu ecuația matriceala z = ABx + Apzq obținem: Czz Concatenam constrângerile pe intrare: CuUi Notând Q = diagțCh, ■■■ , QN) si R = diag(R0, ■■■ , Rn-\), atunci funcția obiectiv devine: = (ÂBx+ApzQ-zref}TQ(ĂBx+ApzQ-zref} + (x-uref}T R(x-uref unde reamintim ca x = [ujf ■ ■ ■ zref = [(z[ef)T ■ ■ ■ (zr^f)T]T si => Forma particulara de programare patratica avand numai constrângeri de inegalitate (matrici dense) microcontroler dsPIC30 (16-bit) Bluetooth senzori infrarosu camera video CMOS (rezoluție 640 x 480) senzor ultrasunete accelerometru 3D • Dezvoltat de EPFL (detalii: ) • Laboratorul de “Optimizare si Control Distribuit” al Depart ACSE are 4 roboti si masa speciala pentru teste (Sala ED205) • Model continuu liniar simplificat al robotului E-Puck (restricționat la deplasarea înainte): • y distanta parcursa, 0 unghiul de viraj, r raza roților, / distanta roata-centru de greutate, ui si U2 vitezele unghiulare ale roților • Notând starea z = [y 0]T si intrarea u = [ui U2]T avem: unde Az = 0 G R2x2 si Bu = 2 21 2 r 21 Discretizam formularea continua prin metoda Euler Pentru pasul de discretizare At avem Zt+1 = (/2 - AtĂz) zt + AtBuut, Alegem At = 0 5 sec si obținem: Problema de control optimal: • Formulare QP rara a problemei de MPC de orizont A/ = 2: min Zi,Uj s l : unde z = [zi Z2], u = [uo ui], Q = diag(Jz, I2), R = diag(Q-ll2, O I/2) • Reformulam următorul QP in variabila x = zf z2r]r: min s L: Metoda de punct interior presupune transformarea echivalentă a problemei QP intr-una fara inegalități: min (1) s L: Ak = b, unde Cj reprezintă linia j a matricei C Metoda de punct interior: se dau un punct inițial x strict fezabil, t > 0, a 0 și parametrul m numărul de inegalități; cât timp mr > e repetă: calculează x(r) ca soluție a problemei (1) pornind din x; actualizează x = x(r) și t = ar 1 5r 1 - 0 5- 0 -0 5- -1 - "150 20 0 8 0 6 Poziția z Referința 40 60 80 100 -viteza unghiulara 1 viteza unghiulara 2 40 60 80 100 Formulam problema de control optimal a menținerii in poziție verticala y Structura: suport de masa pendul format din bila de masa m; tija de lungime /; Vectorul de stare zel4 se compune din: unghiul tijei cu verticala, viteza unghiulara, poziția pendulului (suportului) pe axa Ox, viteza sa Modelul liniar discret al sistemului unde ut G R reprezintă corecția deplasării orizontale, ' 1 0259 0 504 0 0 ' 0 0013’ 1 0389 1 0259 0 0 -0 0504 — -0 0006 0 1 0 05 Bu = 0 0006 —0 0247 -0 0006 0 1 0 025 : menținerea tijei suficient de aproape de verticala (in orice moment avem —10° 0, formulam problema de control optimal pe un orizont finit A/: 52 2z'TQoz/ + 52 2Rou^ " ' i=l i=0 s L: z0 = z, z/+i = AzZj + BuUi @min — (z/)l — @max VI = 0, , A/ — 1, unde z este starea inițiala a pendulului, iar matricele din cost 1 0 0 0' 0 0 01 0 0 0 0 10 0 0 0 0 01 si /?0 = 10 Formularea (QP) rara fara eliminarea stărilor: min -xtQx x 2 s L: kx = b Cx Observam ca la pasul t = 2 constrângerea pe unghi este activa si după pasul t = 15 unghiul tijei cu axa verticala devine 0 = 0 Total number of Websites 800,000,000 Websites 600,000,000 400,000,000 200,000,000 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 An vs Număr pagini web http: //www internetlivestats com / Doar in 2013, numărul paginilor web a crescut cu 30%; Problema centrala a motoarelor de cautare: selectia/ierarhizarea ("ranking”-ul) surselor de informație funcție de relevanta raportata la obiectul cautat PageRank: clasificarea unui număr uriaș de pagini web; Rețeaua paginilor se reprezintă prin intermediul unui , unde nod = pagina, iar muchie = link Ponderea pjj = probabilitatea ca la o navigare aleatorie sa se ajunga din pagina / in pagina J Matrice de adiacenta E e Rnxn: Ejj = ptj > 0 daca intre / si j exista muchie sau Ejj = 0 daca / si J nu sunt legate Matricea E este stocastica pe coloane (i e Ejj = 1 V7), deci valoarea proprie maxima in modul este 1; Determinați vectorul propriu x e astfel incat Ex = x, eTx = 1, x > 0 (x vectorul probabilităților de accesare pagina /) Formulare ca problema de optimizare QP: min -IlEx — x||2 xeR" 2 s L: eTx = 1, x > 0, unde e = r Formulare echivalenta neconstransa via penalitate (alegand t > 0 suficient de mare): 1 T min -IlEx — xll2 H—(eTx — l)2 xeR" 2 II 2V ’ Pe baza teoremei Peron-Frobenius, asiguram constrângerile de inegalitate x > 0 in cazul neconstrans Uzual , restrânge aplicabilitatea algoritmilor la cei de ordinul I ( datorita complexității scăzute per iterație) Aplicam pentru rezolvarea ambelor formulări (/') Formularea constrânsa: min -IlEx — xll2 xeR" 2 s L: eTx = 1, x > 0, Metoda de Gradient Proiectat: xk+1 = [x/ 0} si a > 0 Convergenta metodei de gradient proiectat (n = IO2) Pentru proiecția pe mulțimea simplex An aplicam algoritmi cu convergenta in (necesita O(n) operații) {ii) Formularea neconstransa: 1 T min -||Ex — x||2 + ~(eTx — l)2 xeR" 2 ii 2V 7 Metoda Gradient: x/j+i = Xk — aVF(xk,r) Iterații (k) x Iteratu (k) Curba de convergenta a metodei de gradient pentru problema Google cu n = IO3, t = 50 (stânga); dependenta convergentei de parametrul r (n = IO3) (dreapta) Tehnicile de (“pattern recognition”) se ocupa cu identificarea clasei din care un obiect studiat face parte etapa de antrenare: se acumulează un set de observații (sau instanțe) cu apartenența la categorii cunoscuta; etapa de clasificare: se dezvolta modele matematice (estimatoare) cu rolul de a clasifica un nou obiect cu apartenența necunoscuta; Aplicații: recunoașterea email-urilor de tip spam sau malware; recunosterea vocii/fetei; detecția de tipare in cadrul unei imagini; recunoașterea scrisului de mana • Dispunem de un set de puncte recunoscute a priori-, {y/}^ (e g puncte de culori diferite) Pentru fiecare punct y,- cunoaștem clasa din care face parte: q cu valoarea +1 daca y,- este de culoare roșie sau ci = —1 daca este de culoare albastra • Principiu antrenare: Se determina un model matematic (hiperplan) astfel incat orice y,- ce satiface ary b + £ este de culoare albastra • Presupunem ca datele {y/}^ sunt liniar separabile si selectam doua hiperplane marginale descrise de: aTy — b = 1 si aTy — b = —1 ce nu conțin puncte intre ele • Problema: Determinați (a, b) ce maximizează distanta dintre cele doua hiperplane, descrisa de 2/||a|| Etapa de antrenare se reduce la problema de optimizare: 1 min - aeR",beR 2 2 s L: unde a si b reprezintă parametrii hiperplanului, iar Ci indica clasa (culoarea) obiectului y,- Problema de optimizare convexa patratica avand numai constrângeri de inegalitate => metodele de optimizare constrânsa pot fi utilizate (e g metoda gradient proiectat - dificil de proiectat pe un set mare de semispatii; metoda de punct interior, etc ) => Principiu clasificare: Hiperplanul cu parametrii (a, b), obtinut in etapa de antrenare, se utilizează pentru identificarea unui nou obiect z; daca atunci obiectul aparține clasei 1, altfel daca aparține clasei 2 • Reprezentarea numerica a imaginilor: fiecare pixel este definit de o valoare (e g intre 0 și 256) ce conține culoarea acestuia • Consideram imagini mono-colore de dim 7 x 7 ale cifrei 7 unde pixelii sunt reprezentat! de nivele de gri cu valori intre 0 si 5 • Problema: Determinați daca intr-o imagine data se afla cifra 7 • Etapa de antrenare: acumulam un set de imagini de antrenare ale cifrei 7 in diferite poziții (clasa I) si imagini aleatorii complet diferite de cifra 7 (clasa II) • Imaginii i se asociaza un vector y; E N49 (cu valori intre 0 și 5) si indexul c; al clasei (dacă Q = 1 atunci conține cifra 7, daca Cj = —1 atunci imaginea este aleatorie) • Determinam hiperplanul optim cu parametrii (a, b) Etapa de clasificare: • Pentru imaginile de test (de mai jos) calculam valoarea hi perpla nu I ui: aTy — b /), i = 1, ■ ■ ■ , m funcții de baza p7(t) = t7_1 pentru j = 1, ■ ■ ■ , n matricea A cu intrările A,j = (matrice Vandermonde) /i ■■■ A= 1 t2 k 1 tm t2m ■■■ t"-1 ) • Exemplu: aproximam funcția f(t) = 1+41tQt2 cu un polinom pe intervalul avand la dispoziție 100 de date Regresie liniara (A G R100xc/) cu polinoame de grad d = 1,2,3,4 avand eroarea ||4x* - b\\ = 0 135, 0 076, 0 025, 0 005 Linie continua g(t), linie punctate Pd(t) pe Consideram un sistem intrare-iesire: Z Pentru sistemul considerat dispunem de 40 de măsurători intrare-iesire y(t)}: Dorim sa aproximam sistemul printr-un model intrare-iesire de forma ARMA: ymodel(t) = Xlu(t)+x2u(t-l)+x3u(t-2)+X4u(t-3)+x5u(t-4) Determinarea setului de parametri x = [xi ■ ■ xs]r ai modelului poate fi realizata prin rezolvarea unei probleme patratice neconstranse: min ||Ax — b||2 cu A G R36x5 xSR5 b=[y(5) y(6) - y(40)]r / "(5) "(6) "(4) "(5) "(3) "(4) u(2) u(3) "(1) "(2) A = u(7) "(6) u(5) u(4) "(3) \ "(40) u(39) u(38) "(37) "(36) } Răspunsul modelului estimat: Răspunsul estimat si răspunsul real Tomografie computerizata = tehnica ce folosește raze X (sau alte tipuri de radiații) pentru a produce imagini 2D/3D ale interiorului obiectului scanat Procedura de funcționare: 1 Se achiziționează o serie de , din diferite unghiuri, ale obiectului scanat; 2 Prin intermediul proiecțiilor obținute, se reconstruiește interiorul obiectului cu ajutorul unui algoritm iterativ; www mathworks com In majoritatea cazurilor, radiațiile folosite sunt daunatoare; de aceea se urmărește achiziționarea unui număr minim de proiecții • Fie x G Rn imaginea interiorului de reconstruit • Pentru reconstrucție, dispunem de diferite măsurători liniare (proiecții) ale imaginii x: b; = A;x i = 1, ■ ■ ■ , m • Notam vectorul proiecțiilor b G Rm si A = [4-f ■ ■ ■ e Rmxn matricea de achiziție • Imaginea interiorului reprezintă soluția sistemului liniar (subdeterminat deoarece sunt mai puține măsurători m decât dimensiunea imaginii n): • Reformulare in termeni de problema CMMP: unde de obicei se alege a = 2 sau a = 0 sau a = 1 Alegem a = 0 V 1 pentru a induce o reprezentare rara a imaginii (vectorul soluție) Se dorește o reprezentare rara a imaginii, deoarece aceasta permite: compresie ușoara; algoritmi rapizi pt procesare; memorie de stocare mica; eliminarea ușoara a zgomotului; • Pentru rezolvare, cu a = 2, se poate utiliza metoda gradient sau metoda gradientilor conjugați! Definim o măsură a raritatii ||x||o =» {i : x,- 0} si ||x||i = ^2,- |x,-| si apoi problemele de optimizare asociata sistemului +x = b: min llxllo s l Ax = b min llxlli s l Ax = b xeR" xeR" relaxare convexa la^a-l • Mutual coherence: u(A) = max ,, îmi u l f diferenția bila xeR" Exemplu: sparse data fitting (aproximarea unui set de date cu funcție liniara): set de date (a,-, b,) a i b,- depinde linear de a,-: b/ = ajx + et \/i = 1 : m cu e,- o eroare datorata procesului de sampling => e g CM MP rar min fT(x) (= ||Ax — b||2 + r||x||i) => A G Rmxn xeR" Aplicații: Magnetic Resonance Imaging (MRI): instrument medical imagistic folosit pt scana anatomia si fiziologia unui organism Image inpainting: technica de reconstrucție imagini degradate Image deblurring: instrument de procesare imagine pt indeparta neclaritatea cauzate de fenomene naturale (mișcare) Genome-Wide Association study (GWA): comparare ADN intre 2 grupuri de persoane (cu/fara o boala), in scopul de a investiga factorii de care depinde o boala min fT(x) (= f(x) + t||x||i) => f diferenția bila xeR" Aproximam f cu o funcție patratica si păstrăm nealterat t||x||i: xk+1=arg rn\nqT(x-xk)(=f(xk)±(^f(xk) x-xk}±h\x-xk\^k±T *lli) Hk = Lln unde L constanta Lipschitz a lui Vf => metoda gradient (accelerat): Nesterov’07 Hk = diag(/_i, ■ ■ ■ , Ln), unde L,- constanta Lipschitz a lui V,-f => metoda gradient pe coordonate: Nesterov’10, Necoara’ll Hk = ^2f(xk) si inlocuieste funcția nediferentiabila ||x||i cu aproximarea diferentiabila (Huber) |x,| o ^/x2 + M2 — /' => metoda Newton: Gondzio’15 Probleme de data fitting cer analiza Big Data (terabytes of data): trilioane de varibile (n) => Cray XC30 MPP supercomputer (ARCHER in Edinburgh este al 25-lea in top 500): 118 080 cores cu performanta 1 642 TFIops/s on LINPACK benchmark ARCHER in Edinburgh (25 in top 500): n = IO12 (dim trilion) Comparație metode: gradient pe coordonate paralel (PCDM), gradient accelerat (FISTA), Newton (pdNCG) cu gradient conjugat pt rezolvare aproximativa sistem liniar si quasi-Newton (PSSgb) Definim o măsură a raritatii ||X||o = rang(X) si ||X||* = 52,-cr,- si problemele de optimizare asociate funcției liniare A : Rmxn —> Rp: min ||X||0 s l A(X) = b => min ||X||* s l X(X) = b X^KmXn X£RmXn relaxare convexa • Restricted Isometry Property de ordin r: cea mai mica constanta 5r{A) a i pt orice matrice X de rang cel mult r avem: 1 - 6r(A) 1 Atunci unica soluție a problemei originale neconvexe este matricea de rang cel mult r satisfacand A(X} = b • Presupunem ca r > 1 satisface d$r = o imaginea data (40% elemente lipsa) si imaginea recuperata Video, sunet sau imagine, depind de mărime si rezoluție, si deci pot ocupa mult spațiu pentru a le stoca sau transmite E g , compresia de imagine se imparte in 2 clase: si compresia lossless: toata informația fișierului original poate fi recuperata după decompresie (formatul GIF este un exemplu), important in documente text unde pierderea de informație duce la imposibilitarea citirii textului compresia lossy: informația redundanta este permanent stearsa (formatul JPEG este un exemplu), important in fișiere unde o pierdere de informație mica este imposibil de detectat de om Exista multe posibilități de compresie a semnalelor, dar avantajul oferit de DVS consta in flexibilitatea sa, deoarece poate acționa pe orice matrice de dimensiune m x n In general, procesele inginerești implica achiziția/ prelucrarea/ comunicația unui număr de semnale/date modificări nedorite ale semnalelor (zgomot) rezultate din: Achiziția semnalelor (e g imagine, video, audio) Comunicații imperfecte Exemplu: semnale cu detalii excesive au variație totala mare (i e integrala gradientului absolut al semnalului este mare) ==> reducerea zgomotului echivalenta cu reducerea variației totale Recuperarea perfecta a semnalului original este Concluzie: căutăm “cea mai buna” aproximare a imaginii coru pte Aproximarea se determina utilizând algoritmi de optimizare (pentru niveluri rezonabile de zgomot avem recuperare “aproape” perfecta) Modelare pentru semnale 2D (tip imagine): Asociem imaginii o matrice X G Rmxn sau un vector x G Rmn Dimensiunile (m, n) reprezintă numărul de pixeli pe linii si coloane Ex : Pentru o imagine cu dimensiuni 640 x 480 pixeli asociem X g R640x480 O imagine se reprezintă numeric printr-o matrice 1 Aproximare matrice cu matrice de rang mai mic via DVS: fie A G Rmxn cu rang(A) = r si A = UY VT = căutăm matricea A avand rangul p n — p Apoi: dim Spân {vi, • • • , vp+i} = p + 1 Cele 2 spatii se intersectează: 3z : ||z|| = 1, Bz = 0, z G Span{vi, ■ ■ ■ , vp+i} p+1 p+1 (A - B)z = Az = ^2 ^iUiV^z = z)uî (z ± vp+2) /=1 i=l p+1 Cum u^uj = 6ij avem 11(A-B)z112 = £ct2(v,rz)211Ui112 > a2p+111z112 i=i M - B\\ = „max ||(A - B)x\\ > ||(A - B)z|| " > ap+1 = ||A - Â|| Fie imaginea originala 359x371 pixeli (Melancolia de Durer, 1514): pe care o putem scrie ca o matrice A de dimensiune 359 x 371 (fiecare intrare (/,j) reprezintă intensitatea de gri a pixelului (/,j) cu valori intre 0 (negru) si 255 (alb)), care poate fi descompusa via DVS ca: A = UYVt, unde U este 359 x 359, E este 359 x 371 si V este 371 x 371 Imaginea originala are 359x371 = 133 189 pixeli Dar matricea A poate fi scrisa ca suma de componente principale: A = (TlUivT + ff2U2^2 + H anunv^, unde fiecare matrice utvT (componenta principala) de rang 1 este de marimea matricei originale A Dar pentru ca valorile singulare cr2 > • • • > crn > 0, este posibila o compresie semnificativa a imaginii atata timp cat spectrul valorilor singulare are doar cateva intrări mari si restul sunt mici! Spectrul lui A conține doar 100 — 200 componente principale max l Putem deci reconstrui fidel imaginea folosind doar o submultime de componente principale De exemplu, putem recupera o imagine asemenatoare cu doar de componente principale folosind comenzile din Matlab: [U,S, V] = svd(A) si B = U(:,l : 20) S(1 : 20,1 : 20) V(:,l : 20) Imaginea B folosește doar din memoria imaginii A\ 20 x 359 + 20 x 371 + 20 = 14 620 vs 359 x 371 = 133 189 pixeli Putem recupera o imagine fidela cu doar de componente principale, in loc de 359 cat are imaginea originala (i e 50% informație din imaginea originala): [t/,5, V] = svd(/\) si B= U(:, 1 : 200)5(1 : 200,1 : 200) V(:,l : 200) Consideram o imagine originala 500 x 800 pixeli de rang 500 dispunerea valorilor singulare - prima figura (scala log), putem observa (figura din dreapta) cata informație este deținuta in primele p valori singulare, facand raportul: in imagine! Deși imaginea originala are 500 x 800 pixeli, putem recupera o imagine fidela cu doar de componente principale (10%) % citește imaginea si transform-o in negru si alb: tiger = rgb2gray(imread(tiger jpg)); % resample imaginea: tiger = im2double(imresize(tiger, 0 5)); % calculează DVS-ul imaginii: [L/,S, V] = svd(tiger); % reprezintă valorile singulare (scala log) si informație %: sigmas = diag(S); plot(log(sigmas)); plot(cumsum(sigmas)/sum(sigmas)); % arata imaginea originala: imshow(tiger); % calculează aproximările de rang mai mic: ranks = ; ns = length(sigmas); for i = l:length(ranks) approxsigmas = sigmas; approxsigmas(ranks(i) : end) = 0; approxs = S; approxs(l : ns, 1 : ns) = diag(approxsjgmas); % calculează matricele de rang mai mic si ploteaza: approxtiger = U approxs V; imshow(appraxt,„er) Deși imaginea originala (faimoasa si frumoasa Lena :)) are 600 x 600 pixeli, putem recupera o imagine fidela cu doar de componente principale (20% din datele imaginii originale):